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Kurzfassung 


Diese Arbeit präsentiert neuartige Simulationstechniken für Spritzgusssimu- 
lationen mit faserverstärkten Polymeren (FRPs). 


Spritzguss ist einer der meistverbreiteten Prozesse zur Massenproduktion von 
diskontinuierlich faserverstärkten Polymerbauteilen. Die Prozessparameter 
(Füllrate, Temperatur, etc.) beeinflussen die Bauteileigenschaften signifikant. 
Für eine adäquate Vorhersage der finalen Bauteileigenschaften muss eine 
Simulation alle Prozessschritte (Formfüllung, Nachdruck, Abkühl- 
/Aushärtungsphase, Abkühlung außerhalb des Werkzeuges) beinhalten. 


Während der Formfüllung hat die Strömungsmodellierung oberste Priorität. 
Das komplexe Matrixverhalten muss unter Beachtung von Scherrate, Tempe- 
ratur und, falls vorhanden, chemischer Reaktion modelliert werden. Die sich 
ausprägende Faserorientierung, die von Strömungsfeld, Faserlänge und 
Volumengehalt abhängt, sollte aus zwei Gründen berechnet werden. Einer ist 
das Ausprägen von anisotropen Material- und somit auch Bauteileigenschaf- 
ten aufgrund der Fasern. Zudem rufen die Fasern auch während der Formfül- 
lung anisotropes Verhalten im flüssigen Material hervor. Auch die Faserlänge 
beeinflusst das mechanische und Fließverhalten des Materials und wird im 
Umkehrschluss durch das Strömungsfeld während der Formfüllung beein- 
flusst. Die Faserlänge hat großen Einfluss auf die Schlagzähigkeit des Bau- 
teils, aber auch auf die effektive Viskosität in Faserrichtung im flüssigen 
Material. Umgekehrt erzeugt das Strömungsfeld aber auch Kräfte auf die 
Fasern, die diese zum Brechen bringen können. Stand der Technik Simulati- 
onen beachten den Einfluss der Faserorientierung und -länge auf das Strö- 
mungsfeld nicht. Diese Arbeit präsentiert einen neuartigen Ansatz, in wel- 
chem Viskosität, Faserorientierung, Faserläinge und Geschwindigkeit 
gekoppelt sind. 


Kurzfassung 


Zur Berücksichtigung der Fasereigenschaften in der Viskositätsmodellierung 
und somit auch in der Geschwindigkeit wird die Viskosität als Tensor vierter 
Stufe, der als Funktion von Matrixviskosität, Faserorientierung, -länge und - 
volumengehalt definiert ist, modelliert. Der Viskositätstensor wird für eine 
homogenisierte Matrix-Faser-Suspension auf Basis von mikromechanischen 
Modellen berechnet. Für die Modellierung des Faserbruchs werden die 
hydrodynamischen Schlepp- und Auftriebskräfte beachtet. Zusätzlich werden 
makroskopische Ansätze zur Berechnung der Faser-Faser Interaktionskräfte 
(Schmier- und Reibkraft) gezeigt und verifiziert. 


Neben der Formfüllung beeinflussen die weiteren Prozessschritte Nachdruck, 
Abkühl-/Aushärtungsphase und Abkühlung außerhalb des Werkzeuges 
ebenfalls die Bauteileigenschaften. Durch das anisotrop visko-elastische 
Verhalten können Verzug und Eigenspannungen aufkommen. Stand der 
Technik Software simuliert diese Phänomene in der Regel anisotrop mit 
linear elastischen Modellen. Diese Arbeit präsentiert einen Ansatz zur Be- 
rechnung von Verzug und Eigenspannungen für FRPs mit duromerer Matrix 
und thermo-visko-elastischen Modellen. Relevante Prozessdaten wie Faser- 
orientierung, Temperatur und Aushärtungsgrad werden übertragen um diese 
in der Verzugssimulation mit zu betrachten. Faser- und Matrixeigenschaften 
werden zur Homogenisierung verwendet und unter Beachtung der Faserori- 
entierung wird ein orthotropes Material definiert. Das Matrixverhalten wird 
als Funktion von Aushärtungsgrad und der Temperatur modelliert. Zusätzlich 
werden thermische und chemische Schwindung beachtet. 


Die vorgestellten Methoden sind für Formfüllsimulationen in der open- 
source, finite Volumen basierten Software OpenFOAM und für die Ver- 
zugsanalyse in die kommerziellen finiten Elemente basierten Software 
Simulia Abaqus implementiert. Numerische Studien verifizieren die Imple- 
mentierung und Methoden. Die Formfüllsimulationen zeigen eine gute 
Übereinstimmung mit experimentellen Ergebnissen, was die neu entwickel- 
ten Ansätze validiert. 
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Abstract 


This work presents novel simulation techniques for injection molding of fiber 
reinforced polymers (FRPs). Injection molding is one of the most applied 
processes for mass production of discontinuous FRP parts. The process 
conditions such as filling rate, temperature, etc. have a significant influence 
on the final part properties. For an adequate prediction of these properties, a 
process simulation has to depict different aspects, including all process steps, 
being mold filling, holding, in-mold solidification and out-of-mold cooling. 


During the mold filling phase, the flow modeling is of major significance. 
The complex matrix behavior must be modeled under consideration of 
shearing, temperature and, if present, chemical reactions. The important 
aspect of fiber orientation, depending on flow field, fiber length and volume 
fraction, should be modeled for two reasons. The first one being, that the 
fiber orientation influences the anisotropic mechanical properties of the 
material, and therefore, the final part’s behavior. Furthermore, the fibers also 
produce an anisotropic flow behavior in the liquid material during mold 
filling. Similar to the orientation, the fiber length influences the flow and 
mechanical behavior of the material and is vice versa influenced by the mold 
filling process. The fiber length is crucial for the mechanical impact strength 
and resistance and the effective viscosity in fiber direction. On the opposite 
the flow field evokes forces on the fibers, leading them to break. In state-of- 
the-art simulation techniques, the influence of fiber orientation and fiber 
length on the flow field is not considered. Therefore, this work presents a 
novel simulation approach where viscosity, fiber orientation, fiber length and 
velocity field are coupled. 


For consideration of the influence of fiber properties on viscosity and hence 
velocity, the viscosity is modeled with a fourth order anisotropic viscosity 
tensor, depending on matrix viscosity, fiber orientation, length and volume 
fraction. The viscosity tensor is calculated for a homogenized matrix fiber 
suspension, based on micro mechanical models. For modeling of fiber break- 
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Abstract 


age, hydrodynamic drag and lift are computed. Furthermore, novel macro- 
scopic approaches for fiber-fiber interaction forces (friction and lubrication) 
and contact points are presented and verified. 


Besides the mold filling, the following process steps holding, solidification 
and out-of-mold cooling also have impact on the final part’s properties and 
geometry. Due to the anisotropic and viscoelastic mechanical behavior, 
warpage may occur, or residual stresses build up. State-of-the-art software 
simulates these phenomena anisotropic with linear elastic mechanical mod- 
els. This work presents an approach to calculate warpage and residual stress- 
es for FRPs with thermoset matrix using thermo-chemo-elastic material 
models. Temperature and curing fields are mapped to be considered in the 
warpage simulation. Fiber and matrix properties as well as fiber orientations 
are used for homogenization to create an orthotropic material model. The 
matrix behavior depends on the degree of curing and temperature. Thermal 
and chemical shrinkage are also considered. 


The presented methods are implemented in the open-source, finite volume 
based software OpenFOAM for mold filling simulations and in the commer- 
cial finite element based software Simulia Abaqus for warpage simulations. 
Numerical studies verify the implementations and methodology. The mold 
filling simulations are validated by comparison to experimental results, 
showing good agreement. 
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1 Introduction 


1.1 Motivation 


Reduction of costs, weight and CO2-emissions are major priorities in almost 
every engineering field. Especially in the mobility sector, systems need to be 
cost-effective and light. Due to this development, new challenges such as 
functional integration and low-density materials are rising. 


The aspect of low-density drives discontinuous fiber reinforces polymers 
(FRPs) in focus of engineering tasks. These material systems offer good 
specific mechanical behavior due to fiber reinforcement in combination with 
low-density of the polymer matrix. Due to combination of fibers, matrix and 
several available additives, the material systems can be specified in regard of 
mechanical and thermal properties, chemical resistance, flame protection, etc. 
Furthermore, FRPs have a high potential for functional integration. The low 
temperature of the liquid material (compared to metals) creates the possibility 
of encapsulation of electric components, metal inserts and other functional 
groups during part manufacturing. 


Manufacturing FRPs in combination with functional integration and mass 
production at low-cost leads to injection molding. Injection molding is one of 
the most frequently applied processes for manufacturing discontinuous FRPs 
with thermoplastic or thermoset matrix. Usually, the process is non- 
isothermal, transient and non-evacuated. Due to these aspects and since the 
process conditions have crucial impact on the final part properties, especially 
within the presence of fibers, a high-quality process simulation supports the 
prediction of the parts properties and geometry. Furthermore, a process 
simulation shows the general manufacturability of a given geometry under 
defined process conditions, considering pressure, temperature, air traps etc. 
Based on process simulations, process conditions can be optimized with 
respect to production time and energy effort. These aspects help to lower the 
development time and costs at an early stage of the engineering task. 


1 Introduction 


Although several approaches are already existing to simulate the injection 
molding process and the flow behavior of FRPs has been in focus of several 
research projects, a state-of-the-art simulation still includes many simplifica- 
tions and neglections. Therefore, this work presents a novel approach for all 
stages of the injection molding process, which takes the influence of fibers 
and in-mold air on the filling process into account and simulates the follow- 
ing process steps anisotropic and viscoelastic. 


1.2 Objectives of the thesis 


The major objective of this thesis is increasing the level of detail for an 
injection molding simulation with FRPs for all process steps. Therefore, two 
points are addressed in mold filling simulations. One is the simulation of the 
in-mold air in addition to the FRP phase. The presence of air can lead to air 
entrapments for complex geometries and has influence on the in-mold pres- 
sure and thermal equilibrium. The second point is the influence of fibers on 
the flow behavior of the FRP and vise versa. 


The co-existing FRP and air needs to be modeled during the mold filling 
simulation. Therefore, the simulation is extended from a single- to a multi- 
phase approach. FRP and air are simulated as immiscible fluids with separate 
flow and thermal behavior. 


By further considering the influence of fibers on the FRP flow behavior and 
enabling an anisotropic flow, the viscosity is extended from a scalar quantity 
to a tensor. The tensor is a function of fiber orientation and length. Since the 
fibers break during mold filling, the fiber length distribution is also consid- 
ered during mold filling and fiber breakage is calculated, depending on flow 
induced forces. 


For a higher degree of detail in simulations of the following process steps 
(after mold filling), an anisotropic, non-isothermal structural analysis is 
performed to predict shrinkage, warpage and residual stresses. Relevant mold 
filling results are transferred and the calculated fiber orientation is used to 
create homogenized and orthotropic materials. 


1.3 Structure of the thesis 


1.3 Structure of the thesis 


1.3.1 Structure of content 


This work is structured along the course of an injection molding process. In 
Section 2 the state of the art is given. Section 2.1 gives a short outline of the 
real process and its different phases and materials. Section 2.2 and Section 
2.3 represent the state of the art for the corresponding simulation approaches. 
Section 1 contains the description of all novel approaches and methods, 
developed during the thesis. Within Section 1 the approaches are verified 
with numerical experiments (Section 4.1) and experimental process data 
(Section 4.2). Section 4.3 gives a short outlook on warpage analysis with the 
novel approaches. Section 5 includes the conclusion and outlook. 


1.3.2 Information about notation 


Within this work, index notation and the Einstein sum convention are used. 
Counting indices (for example coordinate directions) are visualized with 
Italian letters (i,j,k,...). Specifications are represented by indices with 
normal letters (air, FRP, ...). For scalars the specification indices are on the 


bottom (i.e. Xspecification), Since they may be potentiated. Specifications of 


us Ati 
vectors and tensors are represented by indices on the top (i.e. x en 


separate them from the counting indices. All tensors of fourth order in this 
work show right hand, left hand and main symmetry, so they can be visual- 
ized in Voigt notation. For Voigt notation the following order is chosen: 
[X11 X22 X33 X23 %13 X12]. Therefore, a fourth order tensor in Voigt 
notation is given by 


Ty T1122 T1133 T1123 T 113 T 112 
T2222 T2233 T2223 T2213 T2212 
T3333 T3323 T3313 T3312 


Tec 
ijkl è 
A T2323 T2313 T2312 
sym. Tj313 T1312 
Ti212 


2 State ofthe art 


2.1 Injection molding process 


2.1.1 Machine and process overview 


Injection molding is a discontinuous process, well suited for mass production 
of complex shaped parts with high tolerances [1]. The major amount of 
research and industrial applications focuses on thermoplastic injection mold- 
ing (TIM), while this work focuses on thermoset injection molding, named 
reactive injection molding (RIM), since thermosets solidify within an irre- 
versible curing reaction (see Section 2.1.2). One manufacturing cycle, also 
named shot, is divided in plasticization, mold filling, holding, solidification 
(cooling or curing) and ejection. The cycle can be completely automated. 


A general build-up of the injection molding machine and a tool is schemati- 
cally illustrated in Figure 2.1. The machine is departed in a plastification and 
clamping unit. The material is feed in granular form through the hopper. 
Plasticization, transport and mixing is realized with a screw, the plasticization 
is supported by a heating system. During plastification, the screw retracts and 
material is accumulated under pressure in a chamber right before the nozzle 
(screw chamber). In order to inject the material into the mold, the screw 
moves forward to displace the material. The movement is controlled with 
hydraulics and usually a constant velocity or a velocity profile is defined. The 
material is injected via sprue (or runner system) into the mold which are part 
of the tempered tool. After mold filling the material solidifies under pressure 
and the part is ejected. During the end of solidification, the plastification unit 
can already melt material for the next shot to keep cycle times low. 


The two major differences between of TIM and RIM, besides the material, 
are the process temperatures and material plastification. A thermoplastic is 
processed and injected at high temperatures (160 °C to 400 °C) into a cooled 
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mold (10 °C to 160 °C), while a thermoset is plasticized and injected at cold 
temperatures (100 °C to 140 °C) into a heated mold (130 °C to 200 °C) [1]. 
Furthermore, thermoplastics have a more complex screw, containing special 
compression, mixing and shearing units, as well as a non-return valve. 
Thermosets are plasticized with a comparatively simple and compressionless 
screw, having no valve and a constant pitch between the flights [1,2]. The 
simple screw may lead to material or pressure losses during processing, but 
enables unscrewing if the material cures within the plastification unit. 


Clamping unit I Plastication unit 


Tool Screw Hopper 


Chamber 
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Mold 


Cooling/Heating 


Figure 2.1: Schematic illustration of an injection molding machine with tool 


2.1.2 Injection molding materials 


This section will only provide a small overview of injection molding materi- 
als. More detailed information about polymer materials, the molecular and 
atomic structure and the resulting properties as well as further applications 
and relevance in industry can be found in [1-5]. 


Injection molding can be performed with polymers of all three main catego- 
ries (thermoplastics, thermosets, elastomers). The RIM process includes all 
reactive materials, which can be thermosets, elastomers and reactive thermo- 
plastics. Since this work is focused on injection molding with thermosets, the 
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usage of this material group should be considered, although all matrix mate- 
rial independent methods presented within this work are also applicable for 
thermoplastic materials. Even though the focus is on thermosets, some 
information about thermoplastics for injection molding should be provided, 
since TIM is the most important process variant for various industry sectors. 


2.1.2.1 Thermoplastics 


Besides the categories amorphous and semi-crystalline, thermoplastics are 
clustered in commodity (or standard), engineering (or technical) and high 
temperature (or high performance) thermoplastics, which also illustrates their 
typical field of usage. Commodity thermoplastics are mostly known from 
packaging sector due to their low price (< 2 €/kg [5]) and weak mechanical 
performance, which is a positive aspect in this case. Due to these properties 
the usage sector is mainly packaging, for example food and medical products. 
The most prominent thermoplastics in this group are PP, PET, PS and PVC 
[4]. 


As indicated by the name, engineering thermoplastics are mostly used in 
engineering applications. Compared to commodity thermoplastics, they show 
better mechanical properties, especially at higher temperatures (up to 150 °C) 
[5]. Due to the application in engineering tasks, the properties are often 
modified to fulfill special mechanical, thermal, electrical but also optical 
requirements. Therefore, additives and fillers like for example glass (or 
carbon) fibers are added to improve stiffness and strength. Typical applica- 
tions are in the automotive and electronic sector, where these materials are 
used for support structures, housings, etc. The majority of parts is made of 
PA (6 and 66), PC and ABS [4]. 


High performance thermoplastics are characterized by good mechanical 
properties, especially at elevated temperatures. In order to improve the 
mechanical behavior, they may be fiber reinforced. Since costs are often a 
minor aspect for applications with high performance thermoplastics, the 
reinforcement is with carbon fibers, to unleash more lightweight potential in 
the most cases. Typical applications are in the aerospace sector and the most 
used materials are PEEK, PEK and PPS. One of the major material systems, 
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illustrating the lightweight potential of this group, is PEEK reinforced with 
30 %-wt. carbon fibers, available from various manufactures. Such a material 
has a tensile strength comparable to steel and aluminum alloys within a 
density of only about 1.4 g/cm? and can be handled by less than 400 °C [6]. 


2.1.2.2 Thermoset naterials 


Opposite to thermoplastics, thermosets solidify irreversibly and cannot be 
remolten. During curing (solidification) the polymer chains cross-link cova- 
lently, hence the atomic structure is represented by a network and not by 
individual chains, compared to thermoplastics [3]. Due to the covalent 
linkage, thermosets offer good mechanical properties and less tendency for 
creep compared to thermoplastics, only having van-der-Waals forces for 
chain coupling and physical cohesion. 


Thermosets are mainly clustered in four groups: phenol formaldehyde or 
phenolic (PF), unsaturated polyester (UPE), epoxy (EP) and cross-linked 
polyurethanes (PU). Thermosets are typically highly filled, with fillers such 
as calciumcarbonat, wood (flour and fibers) or glass (spheres and fibers). Due 
to the high stiffness and strength PFs, UPEs and EPs are mostly used in 
engineering applications such as automotive parts. In the last few years, PFs 
have also been used for engine parts, caused by the good mechanical proper- 
ties at high temperatures in combination with low cost production [7,8]. 
Epoxies show a low viscosity during manufacturing as well as good adhesion 
and dielectrical properties. Therefore, they are often used for electronic 
encapsulation and as base material for industrial glue. Typical PU applica- 
tions are housings of computers, televisions and other commercial electronic 
products, due to the low cost and high impact strength of this thermoset 
group [1]. 


2.1.2.3 Flow behavior of thermoset materials during processing 


The most important material aspect for mold filling behavior and mold filling 
simulations is the viscosity, defining the materials resistance to flowing. The 
viscosity of a polymer is influenced by temperature, shear rate and degree of 
curing/crystallization. Other aspects such as polymer chain architecture 
(length, side groups, etc.), fibers, pollutions and humidity also influence the 
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material behavior, but are not explicitly mentioned for description of the 
viscosity. 


A general course of polymer’s viscosity for different temperatures and shear 
rates is shown in Figure 2.2. The viscosity is lower for higher temperatures, 
which is a typical behavior of liquids. Polymers often show a Newtonian 
behavior for low shear rates, meaning the viscosity is constant (see Section 
2.2.2.3). For higher shear rates the materials act shear thinning and the 
viscosity decreases. The general viscosity level, presence and width of the 
Newtonian area as well as the change rate in the shear-thinning area depend 
on material and fillers [3]. 


Newtonian , shear-thinning 


Viscosity 


Shear rate 


Figure 2.2: Qualitive illustration of polymer viscosity over shear rate for two different tempera- 
tures 


For processing of thermoset (or elastomer) polymers, the irreversible curing 
reaction needs to be considered for viscosity. Figure 2.3 illustrates the 
viscosity of a thermoset during processing. Since the material is heating up 
during mold filling (see 2.1.1) the viscosity would descent continuously 
during processing, but the cross-linking of the polymer chains has crucial 
impact on chain movement and hence rises the viscosity. The superimposi- 
tion of these two effects leads to a parabolic related course of viscosity during 
processing, with a specific time window, where the viscosity is low and the 
material can be processed easily and energy efficient [5]. 
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Figure 2.3: Qualitive illustration of thermoset viscosity during processing for temperature only 
(red), curing only (blue) and combined (green) 


The complex viscosity behavior of thermoset materials leads to a complex 
flow behavior and mold filling. The flow behavior is shown schematically in 
Figure 2.4. 
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(low viscosity resin) (high viscosity resin) wall contact front 


Figure 2.4: Fiber reinforced thermoset's flow behavior during mold filling. Hot surface layer 
(red) and cold core region (blue) in tool (grey) with corresponding fiber 
orientation, temperature and velocity profile 


Due to the temperature difference between tool and material, a hot surface 
layer and cold core region build up. Since the temperature in the surface layer 
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is high, viscosity of the resin is low. For the surface layer a shear dominated 
laminar flow is observed, while the flow in the core region is more plug-like, 
dominated by elongation stresses [9,10]. The overlay of these flow phenome- 
na lead to a running ahead of the core region, resulting in a region with 
partial wall contact and an unstructured flow front [9-12]. Newer studies also 
show wall-slip effects, which influence the velocity profile, being not zero on 
the wall and hence also influencing viscous stress, fiber orientation, etc. [11]. 
This wall-slip effect is often neglected within simulations (see Section 
2.2.1.5). Fibers align along the resulting forces of the flow field. Therefore, 
the alignment is in flow direction in the surface layer and perpendicular in the 
elongation dominated region [9]. Nevertheless, the fibers influence the 
viscosity and flow behavior (see Section 3.3). Fiber alignment in flow direc- 
tion increases the viscosity in this direction and is therefore contrary to the 
effects mentioned above. 


Some studies on the flow behavior of thermosets also show a qualitive 
distribution of the viscous stress over the different layers [9]. In these cases, 
the viscous stress is zero in the center of the core region (plane of symmetry) 
and increases linear towards the tool walls. The linear increase is only true 
for a Newtonian fluid in a perfect parabolic velocity field (Poiseuille flow), 
which are both not the case for RIM processes and should therefore not be 
assumed. 


2.1.3 Mold filling 


Even though mold filling takes up only a small amount of time in a complete 
injection molding cycle, it is an important part defining various aspects of the 
final material state and part. The mold filling has crucial impact on fiber 
orientation, degree of curing, temperature distribution, air traps and energy 
requirement. Furthermore, insufficient mold filling may degrade the final part 
or damage the machine. 


The material enters the cavity at the inlet and takes up the void volume until 
filled, as shown in Figure 2.5. By default, the cavity is not evacuated and the 
FRP displaces the in-cavity air, which leaves via venting slots, being too 
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small for the liquid FRP to enter. In most cases, the mold filling is volume 
controlled (constant volume flow or profile), when the cavity is nearly filled, 
there is a switchover to a pressure-controlled filling. The switchover can be 
defined differently. Most common methods are position of the screw, hydrau- 
lic or screw chamber pressure. Another likely used criterion is reaching a 
defined mold pressure, dedicated by one or more sensors, ideally placed near 
the end of the flow path. Of course, these criteria can be combined and some 
aspect such as maximum screw position and maximum hydraulic pressure are 
generally present to secure the machine. 


start of filling filling 


Figure 2.5: Schematic illustration of the mold filling process in injection molding. Material 
(green) enters at the inlet and fills the cavity (white) 


The mold may contain functional parts such as metal inserts, threads or 
electric chips, sensors and connectors. These components get fixed or encap- 
sulated during mold filling to functionalize the part. The filled cavity corre- 
sponds to the final geometry of the part, besides eventual deviations due to 
shrinkage and warpage. The following process steps after mold filling are 
described in Section 2.1.4. 
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2.1.4 Process steps following mold filling 


2.1.4.1 Holding stage 


After the mold is filled, the next process step is the holding stage, when a 
high pressure (up to 200 MPa) acts from the plastification unit via the inlet on 
the material. The primary function of the holding pressure is to press out 
retaining air and counteract the shrinkage. In general, the shrinkage decreases 
with increase of holding pressure [13]. The amount of holding pressure 
depends on the material, part geometry and maximum clamping force. 
Similar to the filling rate, the holding pressure may be constant or assume a 
defined profile. The holding time depends on material, geometry and temper- 
ature. The holding stage is performed until the inlet is solidified, after this 
point, the pressure is not present in the cavity anymore. The mold should be 
designed, so sufficient holding stage is performable, depending especially on 
wall thickness and position of the inlet [1]. Amount and time can be estimat- 
ed with simplified equations or determined and optimized with process 
simulations. 


If the glass transition temperature (increasing with the degree of cure) reaches 
the material temperature, residual stresses may build up. Inhomogeneous 
temperature distribution and shrinkage as well as anisotropy due to fibers 
influence the evolution of residual stresses [14]. Residual stresses result in 
warpage after demolding. Several studies investigate these phenomena and 
their relations. Experimental works show the influence of holding pressure, 
material temperature and mold temperature [13,15]. The experimental inves- 
tigations show, that warpage and shrinkage may be decreased by higher 
holding pressure and longer cooling or curing times, but this of course 
increases the energy requirements and process time. Therefore, an optimum 
must be found, which is also addressed by different investigations. Simula- 
tion models offer the possibility to perform a DoE with little effort, as for 
example shown by Huang and Tai [16] and Choi and Im [17]. Other approa- 
ches suggest relations with Kriging surrogate models [18] and artificial 
neural networks [19]. 
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2.1.4.2 In-mold solidification 


The solidification is no serial process step, since the material starts to solidify 
directly after thermal activation or mixing, including also the steps mold 
filling and holding. Nevertheless, a general process cycle contains a final 
solidification step, having the single aim to generate a form stable geometry 
and part. This step is necessary since the inlet is usually not the area with the 
greatest thickness, hence there is still liquid (or high viscous) material left in 
the cavity after the holding step. This step is named cooling for thermoplas- 
tics and curing for thermosets and elastomers. For RIM processes the curing 
step has a crucial impact on the final degree of curing and therefore on the 
final part’s mechanical properties, shrinkage and residual stresses. 


2.1.4.3 Ejection and out-of-mold cooling 


After solidification the part is form stable and is ejected from the tool. The 
ejection is supported by multiple cylinders, pressing the part out. Of course, 
the mold must be designed the way, that the part can be ejected undamaged. 
Position, number and velocity of the ejection cylinders influence the final 
warpage [14]. 


After ejection the part cools down to room temperature and may post cure in 
this period. Since the part is no longer bounded to the tool geometry the 
residual stresses result in warpage. The amount of warpage depends on 
geometry, fiber orientation and volume fraction, temperature distribution, 
degree of cure and residual stresses, resulting from prior process steps. 
Additionally, the cooling in the atmosphere after ejection influences the 
warpage by temperature and thermal shrinkage. To reduce the effect of 
warpage or residual stresses, parts may be clamped and/or tempered after 
ejection. 
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2.2 Mold filling simulation 


2.2.1 Process modeling 


2.2.1.1 Governing equations 


In general, the injection molding process is a transient problem of fluid 
dynamic type. These can be solved analytically in some cases or under 
consideration of simplifications. Nowadays, fluid dynamic problems are 
solved in a computational fluid dynamics (CFD) simulations performed on a 
discretized mesh [1]. Solving a CFD problem for the mold filling simulations 
involves solving of the Navier-Stokes equations postulated in the first half of 
the 19% century [20], including the mass, momentum and thermal energy 
equation. The mass balance is shown in the form 


dp GPs: d(pU;) 
dt Ox; 


=0, (2.1) 


where p is the density, t the time, x; the coordinate vector and U; is the 
velocity vector. The momentum equation is re-presented by 


EUs) MU) i) _ SER Oe 


(2.2) 
dt Ox; Ox; Ox; 


with p being the pressure and 7;, the viscous stress tensor. Body forces such 
as gravity would be formulated with terms like pg;, but are often neglected, 
as in the present work, since they only show minor significance. The thermal 
energy equation is described with Fourier’s Law by 


d(pc,T) 5 d(pc,U;T) 
dt Ox; (2.3) 


OT 
= Ath T + TijDij, 


where T is the temperature, Ap and cp are thermal conductivity and specific 
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heat capacity. The term 7;;D;, calculates the viscous heat/dissipation with the 


strain-rate tensor D;;, being the symmetrical part of the velocity gradient 


j? 
[1,20]. Other terms such as thermal radiation and source terms due to chemi- 


cal reactions are neglected within this work. 


For an isotropic fluid (scalar viscosity) the viscous stress tensor is given by 


alpea a o) 24) 
Ti; = 2n | Di — 55 ôi ), . 
ij i 3ax, 4 

where 77 is the dynamic viscosity and ô;j is the unity tensor. For anisotropic 
fluids, the viscosity is described by a tensor of higher order. Similar to 


anisotropic mechanical behavior, the viscous stress for an anisotropic case 
can be written as 


Tij = NijkıDkı- (2.5) 


The symmetry rules of continuum mechanics also apply for fourth order 
viscosity tensor Nijgı, therefore it contains up to 21 independent entries. 
Nevertheless, at state of the art the anisotropy of FRPs is usually neglected 
and the viscous stress is calculated with an isotropic viscosity as described in 
Eq. (2.4). 


2.2.1.2 Analytical approaches 


Since fluid dynamic problems exist far longer than the possibility of solving 
numerical methods on computers, analytical approaches and solutions exist 
for special cases. Nevertheless, the solution of the Navier-Stokes equations in 
a 3D-case is part of the Millennium Problems, therefore all solutions come 
along with simplification, assumptions and restrictions. 


One of the mentioned special cases is the so called Poiseuille flow, describ- 
ing the laminar and stationary flow of a Newtonian fluid in an infinitly long 
circular tube with a pressure difference between the two ends and radius R. 
The solution was developed by Hagen and Poiseuille in the mid of the 19% 
century [20]. By assuming the velocity to be zero at the walls and the pres- 
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sure to be constant over the tube’s cross section, the flow field is described 
by a parabolic velocity profile as shown in Figure 2.6. The velocity is given 
by 


A 
U(r) = ma — x3), (2.6) 


where Ap is the pressure drop within the distance I. 


Vs 


Figure 2.6: Poiseuille flow velocity profile (green) in a circular tube 


The equation can be modified for non-circular cross-sections [20]. Today the 
Poiseuille flow is still an often-used model for flow description and verifica- 
tion and therefore is used in various scientific publications. 


Another milestone in flow field modeling is the description of a fluid in a 
narrow gap between two parallel plates, presented by Hele-Shaw at the end of 
the 19% century [1,20]. Considering an isotropic, incompressible fluid with 
no body forces and no wall slip the flow field is described by 


0 (s Piy 0 (s oP’) = 0 
Ox, HS Oy, 0%, HS Ox, u 


X3,max 2 , (2.7) 
x3 
Sus = J 7 dx3 


with the fluidity Sys. Here, the narrow gap between the two plates is in the 


17 


2 State of the art 


range 0 < x3 < X3 max, the two plates are at x3 = 0, x3 = X3 max and perpen- 
dicular to the x-direction. U3 is zero in the complete domain. 


It should be noted, that the viscosity is within the integral since it is not 
assumed to be constant and may depend on x3, temperature or shear-rate. For 
a Newtonian fluid, the integral can be solved independent of the viscosity and 
the velocity is given by 


= 2n 0x (%3,max — %3) 


1 dp i 


2. = 2700, (X3,max — X3) 


(2.8) 


resulting in a parabolic velocity profile. 


Due to the thin wall character of many injection molding parts, the solution 
of Hele-Shaw is well suited for injection molding simulation, even though it 
was developed long before the first injection molding processes. Therefore, 
many of today’s 2D and 2.5D mold filling simulation rely on Hele-Shaw. 
Nevertheless, it should be handled carefully, since it is only valid for thin 
areas far away from the inlet and side walls. [1] 


2.2.1.3 Numerical methods 


Today it is standard to solve the governing equations for a CFD-simulation 
under consideration of numerical methods. The most prominent numerical 
methods are Finite Element Method (FEM), Finite Volume Method (FVM) 
and Finite Difference Method (FDM). The principle of these three methods is 
identical, the geometry is discretized into a finite number of elements as 
shown in Figure 2.7. The problem is then solved for the single elements (or 
cells) and combined to a global solution. The simulation mesh can be realized 
with different levels of complexity as illustrated in Figure 2.7. The simplest 
mesh is a 2D-midplane mesh, representing the geometry by one or more 
meshed planes. Of course, also a 1D-approximation is possible, but most 
injection molding simulation software perform at least a 2D-simulation. The 
2.5D-mesh discretizes only the geometry’s surfaces and is also known as 
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Boundary Element Method. Additionally, to the boundary surfaces, a finite 
number of ‘inner planes’ can be regarded as indicated by the blue elements in 
Figure 2.7. 


The most complex variant is the 3D-mesh, representing the complete geome- 
try with 3D-elements. While the FEM is solved on nodes at the corners and 
edges of the elements, the FVM focuses on the cell centers and surface 
fluxes. It must be distinguished between boundary nodes/faces, where 
boundary conditions are applied and internal nodes/cells, where algorithms 
solve the problem. 


Original Geometry 2D-midplane mesh 


2.5D-boundary mesh 3D-mesh 


Figure 2.7: Original Geometry and different possibilities to discretize the geometry in a finite 
element mesh 
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Every numerical method has its individual advantages and disadvantages. 
While FEM is well suited for Langrage approaches, with the nodes being 
able to move, FVM is flux based and therefore well suited for Euler ap- 
proaches and hence CFD simulations. However, most commercial injection 
molding software is based on FEM, due to the numerical efficiency, enabling 
faster solving [21-23]. Nevertheless, injection molding simulation is a CFD 
problem, which is most likely solved with FVM in other cases. Therefore, 
some scientific publications offer FVM-based solutions for injection molding 
simulations [24-26]. Hence, FVM-based simulation approaches are applied 
within this work. 


Besides the three mentioned discretization-based methods CFD problems are 
also solved with Smooth Particle Hydrodynamics (SPH). This mesh-less 
approach splits the simulation domain into single particles considering 
neighbor particles in a defined area for solving. SPH has been successfully 
used to describe FRP flow processes on microscopic scale in [27]. 


2.2.1.4 Influence of mesh size 


The mesh size has a crucial impact on the solution quality and numerical 
effort. One the one side a finer mesh, meaning smaller elements or more 
discrete points, predicts a more realistic solution, but on the other side this 
directly correlates to more numerical effort, since more equations must be 
solved. Nevertheless, deviations from the analytical solution may occur if the 
mesh is not fine enough. A typical case is the approximation of a parabolic 
velocity profile, as it may occur in injection molding (see Section 2.2.1.2). 
The analytical solution and approximations with increasing number of nodes 
are shown in Figure 2.8. 


The parabolic distribution is discretized by three, five and nine nodes with 
linear interpolation. While the three node solution creates great deviations 
and is not able to display the analytical solution, the five node solution 
already shows a quite good agreement. The nine node approach shows only 
small deviations and displays the analytical solution with good agreement. Of 
course, the interpolation between the nodes could also be of higher degree. 
Therefore, the parabolic distribution can be displayed exact by three nodes 
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with interpolation of second order, but the higher order comes along with 
greater numerical effort. In summary, the mesh fineness is always a compro- 
mise between level of detail and numerical effort. 


Since the quality of result depends on the mesh, the development of novel 
methods should include a mesh study, showing a convergence of results for 
different mesh refinements. A convergence can be reached if the approximat- 
ed distribution can be represented correctly (for example Figure 2.8 with 
quadratic interpolation), or if a finer mesh shows no more improvement for 
the results. A multiphase simulation with immiscible single phases (FRP and 
air in this case) always depends on the mesh, since the theoretical interface is 
reconstructed, depending on mesh dimensions [28]. Nevertheless, the flow 
front shows same tendencies and characteristics for adequatly fine meshes, 
although the results are not identical. This aspect is further investigated in 
[26] and Section 4.1.4. 
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Figure 2.8: Approximation of parabolic distribution with different number of points and linear 
interpolation between the points 
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2.2.1.5 Simulation model 


The mold filling simulation is a transient simulation, where the FRP enters 
the simulation mesh at a defined node for FEM or a defined boundary surface 
(inlet) in FVM. While FEM simulations with commercial software are single- 
phase, considering only the FRP, the FVM-based simulations presented in 
[24-26] perform a multi-phase simulation, considering FRP and air. 


The material in-flow is defined as a boundary condition and can be chosen 
similar to the real process, as described in Section 2.1.3. This also applies for 
the switchover point. The regarded geometry usually contains only the cavity, 
not the tool or plastification unit. Figure 2.9 shows different states of filling 
for a multiphase mold filling simulation of an electrical engine housing. 
Although there is a wall-slip in RIM processes, as described in Section 
2.1.2.3, this effect is often neglected within a simulation, as it is within this 
work. A no-slip boundary condition is applied at the walls, justified by the 
need of less parameters and still creating good results [25,26]. 


inlet 
; I 


Figure 2.9: Different states of a multiphase mold filling simulation of an electrical engine 
housing. Cavity filled with air in transparent grey and FRP in red 


The FRP enters the model at the inlet on the top and the simulation is per- 
formed until the cavity is completely filled. One reason for performing this 
simulation is to verify if the part is manufacturable by the chosen process 
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parameters with respect to maximum machine pressure, clamping force, air 
traps, material temperature, etc. Another important aspect is information 
about the material state needed for further analysis, including temperature 
and fiber orientation distribution as well as solidification state. 


2.2.2 Material modeling 


2.2.2.1 Thermal properties 


The relevant thermal properties are specific heat capacity and thermal con- 
ductivity (Eq. (2.3)). Additionally, the reaction enthalpy may be a material 
property, relevant for thermal modeling, but the influence of reaction heat is 
often neglected, as within this work [14,25,26]. This assumption is justified 
by the thin wall character of the molded parts, so the influence of the reaction 
heat is neglectable, compared to the thermal influence of the massive and 
tempered tool. 


The specific heat capacity c, is a scalar property. Independent of polymer 
type, Cp increases with temperature and shows a jump for the phase change 
[29,30]. Kamal and Ryan [29] present measurements of the specific heat 
depending on temperature and degree of cure showing also a significant 
increase of cp during curing. Since only the liquid state is relevant for mold 
filling and the degree of cure has only marginal changes in this process 
phase, it is often assumed to be constant or only function of temperature. 


The thermal conductivity Ap also increases with temperature and jumps at 
the phase change. Also similar to cp, an increase of Ap is shown in [29]. For 
the same reasons compared to specific heat capacity, the thermal conductivity 
is often assumed to be constant or only function of temperature. The thermal 
conductivity may be a multi-dimensional tensor for anisotropic materials. 
The homogenized thermal conductivity can be calculated with the same 
assumption about serial and parallel connections as performed for mechanical 
properties [31,32]. Nevertheless, the anisotropy of An is also often neglected 
during form filling, since the values are generally low (< 0.5 W/(m:K)) for 
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most polymer composites and the mold filling is performed within a short 
period of time. 


2.2.2.2 Curings kinetics 


The curing reaction, meaning the cross-linking of molecular chains with 
covalent bonds, of thermoset materials has crucial impact on the viscosity 
and hence on the mold filling behavior, even though the fraction of cross- 
linked material is small during mold filling. Therefore, the curing kinetics 
should be modeled accurately for a mold filling simulation of a RIM process. 


The curing of thermosets is an irreversible and exothermal process, which 
can be modeled by mechanistic or empirical models. Mechanistic models 
describe the chemical processes during cross-linking and therefore depend on 
detailed, material specific information, which are hard to determine [14]. 
Empirical models focus on simpler description with less parameters. Here, 
the reaction process is described by a simple differential equation, depending 
on a few empirical parameters, degree of cure and in some cases temperature. 
The simplest possibility to describe the reaction is given by 


dc 
aoe —c)", (2.9) 


with c being the degree of cure and n the reaction order in this case. One 
disadvantage of this approach is, that the maximum change of cure is always 
at c = 0, which is no typical behavior of industrially used thermoset materi- 
als [14]. 


Kamal and Sourour [33] extended this approach by 

dc 

—~c™(1-—c)", 2.10 

ne (1-c) (2.10) 
with an empirical parameter m in this case, so the change rate is zero at the 
beginning of the reaction. Finally, the curing kinetics are calculated with the 


so called Kamal-Malkin (or also Kamal-Sourour) kinetics model [33], deter- 
mining the rate of change as 
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dc 

Fr (K, + Kyc™KM)(1 — c)"KM, (2.11) 
with the empirical parameters Myy and Nyy. The parameters Kı and K, are 
often given by an Arrhenius type approximation in the form 


—Exmi, 
Ky 2 = Axmı,2 ' EXP (Fer) (2.12) 
g 


with empirical factors Agmı and Axm>, activation energies Exy, and Ekm2 
and universal gas constant R, [14]. 


The literature offers more complex models for better descriptions, if the 
temperature is for example near the glass transition temperature [14,34]. 
Nevertheless, the Kamal-Malkin model is often used and well established for 
curing kinetics model for different materials and processes [23,25,26,34-36]. 
Within this work the Kamal-Malkin kinetics model (Eq. (2.11)) in combina- 
tion with Eq. (2.12) is used to model curing kinetics during the complete 
process simulation. 


2.2.2.3 Viscosity properties 


The viscosity must be modeled to determine the viscous stress tensor, which 
is needed to solve the momentum equation as shown in Eq. (2.2). The sim- 
plest approach to model the viscosity is to assume the FRP as Newtonian, 
resulting in a constant viscosity 


n = const. (2.13) 


Meyer et al. [37] and Sommer et al. [38] present good results for simulating 
compression molding with Newtonian matrix behavior. Nevertheless, poly- 
mers show a clearly non-Newtonian behavior, with a viscosity depending on 
temperature and shear rate, described in2.1.2.3. 


The simplest description of a shear rate depending viscosity in case of simple 
shear load is given by the power-law approach [39], determining the shear 
stress T by 
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T = Kp - y™Pt, (2.14) 
and with t = ny the dynamic viscosity by 
n = Kp, : ýD, (2.15) 


with y being the scalar shear rate, defined as ý = ?/D; jij. Kpu and np, are 
constant, positive model parameters in this case. The power-law approach 
models shear thinning (n < 1), shear thickening (n > 1) and Newtonian 
behavior (n = 1). One disadvantage, is the prognostication of a zero-shear 
stress and zero viscosity for zero-shear rate, which may be true for the shear 
stress, but not for viscosity. To overcome this, the approach can be extended 
to describe a so called Herschel-Bulkley fluid [40] with the viscosity given 
by 

ee DE, 2.16) 

Toy +K-y Ý 2Yo 


where Yo and Ty are material specific parameters. For T < Tọ the material 
behaves as a solid. Y, indicates the shear rate for a change from Newtonian to 
shear thinning behavior (cf. Figure 2.2). Today’s most common approaches 
for viscosity modeling of polymer melts are of Cross type. Cross [41] de- 
scribes the viscosity by 


= L 
1-ncwLF ’ 2.17 


with zero-shear viscosity 79 and T* indicating the critical shear stress for 
beginning of shear thinning behavior. Ncw, is a material specific parameter. 
In order to respect the temperature dependence, No is often described by the 
Williams—Landel—Ferry (WLF) equation like 
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2 Acwırı(T — Tg) ) 
Acw r2 + (T a To) (2.18) 
T, = Dewur2 + DewuraP 


No = Domme ( 


Acwır2 = Acwir3 + DewırsP 


where Tg is the glass transition temperature and 


Acwırv, Acwır3 Dewırı, Dewır2 and Dew ir3 are material specific fitting 
parameters. For most materials, it is assumed to be Dew ir3 = 0, so T, and 
Acwır2 are constant. The combination of Eq. (2.17) and Eq. (2.18) is the so 
called Cross-WLF viscosity model, being the most used viscosity model for 
injection molding simulations with thermoplastic materials. 


Although the Cross-WLF model is well suited for thermoplastic materials, 
there are better approaches for thermoset materials. Castro and Macosko [42] 
reported a better fitting by describing the temperature dependence with an 
Arrhenius equation with an additionally term to take the curing kinetics into 
account. For thermoset material the today’s most common viscosity model is 
the so-called Castro-Macosko (CM) model, combining a Cross type equation 
with Arrhenius equation and curing dependence by 


u No ( Cg o 
= „n 1-NcM = 
NoY CoC 
1+ (127) 5 f (2.19) 
Tem 
No = Bemexp (=) 


with material specific fitting parameters Bem, Cemi» Ccmz: Necm: Tem and Cg 


being the gelation conversation point. 


2.2.2.4 PvT modeling 


The aim of pressure-volume-temperature (PvT) modeling is the relation 
between pressure, specific volume (or density) and temperature, also often 
titled as equation of state. The simplest way to describe this relation is to 
assume incompressibility, so 
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p = const. (2.20) 


or equivalent €,, + £22 + €33 = 0, with &;; being the strain rate tensor. This 


assumption is valid for polymer injection molding and often chosen in 
different studies [24,37,38,43]. 


Since this study includes the air within the process simulation, this phase 
would also be considered incompressible, which is an unsatisfying simplifi- 
cation. The PvT behavior of gases is often given by the perfect gas law so 


p 
TR? 


p= (2.21) 
with R, as specific gas constant. The best related description of liquids in this 
case is to assume a perfect fluid, determining the density with 


p 


p= 


Here, po is the density at T = 0 K. 


During the history of polymer processing several PVT models with different 
complexity and experimental effort have been created. An overview of 
different models in relation to different polymers is given by Rodgers [44] 
and Junior et al. [45]. Both studies name the so-called Tait model (or Tait 
equation) as most used and best fitting model for a variety of polymers and 
pressure/temperature scales. The Tait equation represents the specific volume 
v by 


v(p, Ta 


= v (T) h- Cln (1 poe (2.23) 


+ ‚T). 
C, is an universal constant with value 0.0894. The Tait model is also called 
two-domain Tait model, since the formulation changes, depending on the 
volumetric transition temperature, defined as 


28 


2.2 Mold filling simulation 


Tı(p) = bs + bep. (2.24) 
The subfunctions, needed for Eq. (2.23) are given in Table 2.1, wherein Tr = 


T — b,. All mentioned b-parameters within this Section are purely empirical 
and material-specific fitting parameters. 


Table 2.1: Definition of subfunctions for Tait model (Eq. (2.23)) 


Function Tp) <T Tp) >T 

vo(T) Dim + b2mTr bis + b2sTr 

B(T) bzmeXp(-b4mTr) b3,exp(—b4sTr) 
v,(p,T) 0 b exp(bsTr — bop) 


Wang et al. [46] further improved the model, by describing the b-parameters 
as function of the temperature change rate and not as constants. 


Due to a lack of available PvT data, the FRPs within this work are described 
nearly incompressible. Since multi-phase approaches are performed, a com- 
plete incompressible simulation would be unsatisfying with respect to the air 
phase. Therefore, the air will be assumed as perfect gas (Eq. (2.21)) and the 
FRP is assumed to be a perfect fluid (Eq. (2.22)), where R, is chosen very 
high, so p/(TR,) = 0 and p = po. 


2.2.3 Fiber orientation modeling 


2.2.3.1 Movement ofa single fiber 


The modeling of fiber orientation has been in the focus of research for several 
decades. Today’s models are based on Jeffery’s work to describe the orienta- 
tion change of a single ellipsoid in a Newtonian fluid with the so called 
Jeffery’s equation [47]. The orientation of a single fiber is described by the 
nominated orientation vector q;, as shown Figure 2.10. The orientation can be 
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formulated in Cartesian coordinates or with the orientation angles p and 0, as 
shown in Figure 2.10. The Jeffery’s equation for q; is 


dq; _ w= 1 


u nude 72 +1 (Dya; -= qi(aiDij;)), (2.25) 


with the vorticity tensor W;;, being the unsymmetrical part of the velocity 


J’ 
gradient and fiber aspect ratio rç, defined as 


re = Le/ de, (2.26) 


where Le and dẹ are fiber length and diameter. 


Figure 2.10: Orientation of a single fiber (blue) represented by vector q; or angles 0 and p 


Fiber matrix suspensions are clustered in three regimes: dilute, semi-dilute 
and concentrated. The regimes are defined by fiber volume fraction ®p or 
nL;, with n being the number density of fibers in the suspension. In dilute 
suspension it is Øp «1 and nls K 1 [48] or Oe K 1/7 [49]. The space 
between the fibers is large and fibers rotate without influence by their neigh- 
bors, Eq. (2.25) is valid. The stationary solution in this case is a periodic 
rotation of the fiber. In semi-dilute suspensions (®; « 1 and nL; > 1 [48] or 
1/7? «K d;.«1/r for isotropic orientation [49]) fiber-fiber interactions 
occur. If also &; is no longer small, a concentrated regime is reached, fiber 
interactions significantly influence the movement and fibers are no longer 
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free to move independently. At state of the art, the movement of fibers in 
semi-dilute and concentrated regimes cannot be described analytically. 


Folgar and Tucker [50] developed a semi-empirical model to take fiber 
interactions into account. The Folgar-Tucker-model is an extension of the 
Jeffery’s equation with an empirical interaction coefficient Cj. The stationary 
solution is a final orientation, depending on Cı. The approach still describes 
the evolution of a single fiber over time. Applying this model, meaning 
calculate the movement of every single fiber within a mold filling simulation, 
creates a huge numerical effort, which is unacceptable at state of the art. 


2.2.3.2 Macroscopic orientation modeling 


Advani and Tucker [51] present a homogenization scheme, describing the 
fiber orientation evolution with orientation tensors to reduce the numerical 
effort. The exact formulations of the second and fourth order tensors A;; and 


A;jkı are 
Aij = | Yla)aiajdq (2.27) 
and 
Ayjkı = | Y(9) 9:9 j9KU94, (2.28) 


where w(q;) (or W(@,@)) is the probability density function for a specific 
orientation. Inclusion of A;; into the Folgar-Tucker model leads to 


Je | Era = —(WiAnj — AirWgj) 
ie Sl 
+1 


EN (2.29) 


(Dir Ar; + Ai Duy — 2DriAkıj) 
+2C,7(6;; — 34;,)- 
This approach builds the base for several actual orientation models with 


different focuses [52-56]. 
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Today’s most commonly used models are the reduced strain closure (RSC) 
model by Wang et. al. [52] for slow orientation due to high fiber volume 
fractions and the RSC-ARD-model (anisotropic rotary diffusion) by Phelps 
and Tucker [53] for additionally long fibers. The models become more 
complex to consider more effects by need of more empirical parameters. To 
keep the empirical parameters low and since small aspect ratios (7 < 100) 
being in the focus of this work, the RSC model is used. Based on the Folgar- 
Tucker model and by considering the eigenvectors and eigenvalues of the 
orientation tensor, the approach describes the evolution of the second order 
orientation tensor of A;; by 


dt OX, 


r-1 


rê +1 
—2[Aije + (1-«) 

(Lijra v2 M; jmnAmnkı)]Pxı} 

+2KCy(6;; — 3A; ). 


= (Wir Ak; z AirWgj) 


+ {DikAkj + AikDkj 


(2.30) 


The fourth order tensors Lijgı and M;;x, are calculated with the eigenvectors 
v; and eigenvalues A of A;j as described in [52]. The so-called strain reduc- 


tion factor « is empirical and material specific, with x € [0,1]. For x = 1 the 
RSC-model reduces to the Folgar-Tucker-model (Eq. (2.29)). 


2.2.3.3 Characteristics of Aj; 


Since the orientation tensor is a complex construct, information to interpret 
the results is required. The side entries reflect the deviation to the coordinate 
axis, hence it is A;; € [—0.5,0.5] for i + j. A fiber has no front or back end, 
therefore, q; = —q; is valid and A,, € sym. The major entries represent the 
amount of fibers orientated in the respective direction. Due to normalization 
it is tr(A;;) = 1, which leads (in combination with Aj; € sym.) to five 
independent entries [48]. 
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One simple and often used interpretation is to see the eigenvectors as indi- 
vidual fibers with an orientation probability of its corresponding eigenvalue. 
Therefore, transversely isotropic materials can be created by only regarding 
the direction with the highest eigenvalue, or the tensor is used for orientation 
averaging. Further information is given in Section 2.3.2.2. 


Figure 2.11 shows specific fiber orientations with corresponding orientation 
tensor. The three cases quasi-isotropic, aligned and 2D random (Figure 
2.11a-c) are often used orientations to verify novel approaches. 


a) b) 


c) d) 


0 0 o0 0 0 o0 
Aj =|0 1/2 | Ay=|0 1/2 | 
0 0 1/2 0 0 1/2 


Figure 2.11: Different fiber orientations with corresponding orientation tensor. Quasi- 
isotropic/3D random (a), unidirectional/aligned (b), planar isotropic/2D ran- 
dom (c) and fabric like/2D aligned (d) 


Comparing Figure 2.11c and Figure 2.11d illuminates one disadvantage of 
the orientation tensor, since it is identical in both cases, but the orientations 
are clearly different and the mechanical behavior of the final parts would 
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show significant deviations. Similar configurations can also occur for 3D 
orientations. 


2.2.3.4 Determination of C; 


The interaction coefficient represents fiber-fiber interactions, thus it depends 
on @, and fiber orientation [57]. In past works, Cy is determined with expe- 
riments [58,59] or mathematically as function of &; and 7; [48,60]. Heinen 
compared the experimental and numerical approaches in her PhD thesis [61]. 
The results show, that a combination of [48,60] such as 


Ci = 


poe — exp(—0.2244; 77) Pp rp S13 (2.31) 


0.0184 exp(—0.71484, rp) Pir > 1.3 


fits best to the experimental work in [58] and [59]. Within this work Cj is 
calculated according to this approach, which is similar to [61]. 


2.2.3.5 Closure approximations for A;jkı 


The change rate of A;; requires the unknown fourth order orientation tensor 
Ajjx:- Similar to the second order tensor, a change rate equation can be 
determined for the fourth order tensor, which would create numerical effort 
and require an unknown orientation tensor of sixth order [51]. Therefore, 
A;jrı is determined with a closure approximation. 


The first solutions for such closure approximations are the linear approach, 
which is exact for an isotropic distribution, the quadratic approach, which is 
exact for aligned fibers, and the hybrid approach, which is a combination of 
linear and quadratic [62]. Today there are several approaches for closure 
approximations to determine Aij, based on invariants or eigenvalues of Aj; 
[63]. Within this work the IBOF5 closure approximation developed by 
Chung and Kwon [64] is used to determine A;jx,, since was found to be a 


good approximation with quite little numerical effort. 
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2.2.4 Fiber breakage modeling 


2.2.4.1 Modeling of decreasing fiber length 


During processing, or during compounding, fibers may get damaged by 
interaction with the wall or other fibers as well as by hydrodynamic load of 
the surrounding matrix. Damage in this case is similar to fiber breakage, 
hence the average fiber length within a compound or part is at least as long as 
in the initial material or shorter. Since the fiber length has crucial impact on 
the mechanical behavior of an FRP, an adequate prediction of the final fiber 
length distribution is crucial for subsequent material modelling. 


One of the first approaches to model fiber length distributions during com- 
pounding is presented by Shon et al. [65], describing the length evolution by 


dl; 
ae —Kshon (Le — Loo), (2.32) 


with kinetic constant Kspon, determined for different process devices and the 
equilibrium fiber length Ls, representing the shortest possible fiber length. 


Inceoglu et al. [66] extended this approach by considering the specific me- 
chanical energy Ilme so 


dl: zB 
dt = —K” - Ilme ' (Le = Lfo). (2.33) 
Here K” represents a change rate constant, determined by data of a batch 


mixer. 


2.2.4.2 Fiber breaking mechanisms 


The approaches of Shon and Inceoglu may describe the shortening of fibers, 
but in an empirical and linear way. The model parameters may be well suited 
for compounding processes, but the variety of process conditions and geome- 
tries in injection molding is too high to describe the complex process of fiber 
breaking with such simple models. 
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In a first step it is crucial to identify the mechanisms, which lead the fibers to 
break. In principal these are identical to other damage mechanisms, meaning 
a specific strain or stress is reached. However, the material behavior of fibers 
and matrix during processing is complex and it is ambiguous which load case 
or cases are present, so it is unclear if the fiber is breaking due to bending, 
stretching or buckling. Several studies identify buckling as the major, but not 
singular, effect for fiber breakage [47,67-70]. Of course, effects such as 
interaction with walls or other fibers as well as imperfection in material and 
fiber geometry may favor or complicate this process. 


Assuming buckling as breaking criterium, it is necessary to verify that the 
stress field of the fiber, induced by the flow field, puts the fiber in compres- 
sion. Therefore, the orientation between the fiber and the flow field must be 
known. The orientation state, where the fiber is under compression, is known 
as Jeffery orbit [47] and defined by 


Figure 2.12 shows the values of D;jqiq; for all possible orientations, repre- 
sented by a unity sphere in case of simple and normed shear flow. According 
to the color bar the Jeffery orbit is given by the blue area. 


After the load case is known, the acting forces and stress distribution must be 
determined, to evaluate if the fiber is buckling. A first approach is given by 
Burgers [71], assuming the forces to act from the fiber end to the center along 
the fiber axis. In a later work Hernandez et al. [68] extended this approach, 
determining analytical and pseudo-analytical solutions of the force integral 
along the axis, showing a non-linear distribution. This aspect is quite im- 
portant for determining the buckling point, meaning the critical force above 
the fiber buckles. Durin et al. [69] performed simulations for fiber loads with 
different orientations, flow fields and aspect ratios, also based on the model 
of Burgers. Within the study it is shown, that for most cases the buckling 
criterion is reached before the breakage criterium (critical stress). The con- 
clusion is, that for brittle fibers, such as glass or carbon, buckling then will 
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lead to breakage and is acceptable and meaningful as singular breakage 
criterium. 


0.5 
4 
€ 
330 0 
& 
1 
i a 
5s, 
-0.5 


Figure 2.12: Visualization of Jeffery orbit on a unity sphere in simple shear flow. Colors 
represent value of [D; jlaa j, blue regions (negative values) represent Jeffery 
orbit 


Defining buckling as the mechanism for fiber breakage leads to the need of 
defining a point at which the fiber buckles. Therefore, the Euler buckling 
modes, derived by Euler in the middle of the 18" century are a common 
model. The critical load of this purely elastic deformation mode can by 
defined by stress, so 


TC 2 
Obu = Er (4) ; (2.35) 


with Ep being the elastic modulus of the fiber, used by Durin et al. [69]. The 
criterion can also be defined for forces as 


T’ Ede 
Fou = gaz (2.36) 
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used by Phelps et al. [70]. Both works identify the first buckling mode as 
critical, leading to the assumption, that fibers break only at one point, so only 
two new fibers are created. The work of Meyer et al. [27] shows, that at least 
in dilute suspensions also other eigenmodes can be reached and fibers may 
break in three parts. 


2.2.4.3 Fiber length distribution models 


To predict an adequate fiber length distribution within the final part, a physi- 
cally based breakage model considering material properties and process 
condition must be defined. It can be seen in literature, that a crucial part of 
fiber breakage takes place within the plastification unit in case of injection 
molding. Nevertheless, this part of the process is often quite similar and the 
conditions change only slightly, so a known state of fiber length distribution 
at the end of plastification can be used as input for the mold filling simula- 
tion. Two actual models for fiber length distribution modeling with breakage 
mechanisms are given by Durin et al. [69] and Phelps et al. [70]. While the 
first one is built to model breakage in a twin screw extruder and uses single 
fiber orientations, the latter one is focused on mold filling in macroscopic 
injection molding simulations. Therefore, only the approach of Phelps et al. 
will be further discussed. 


In a first step, some assumptions must be defined. The most important are, 
that there is a minimum fiber length, where no more fiber breakage is possi- 
ble and all possible fiber lengths are a multiple of this length. Furthermore, 
the flow field is assumed to be pure shear [70]. Additionally, the force acting 
on the fiber must be determined, to compute if the fiber buckles. Therefore, 
Phelps et al. use the slender-body analysis presented by Dinh and Amstrong 
[72], representing the acting force by 


Inulr 
Fe= = (-D;;q:q;), (2.37) 


where ¢ is a dimensionless drag coefficient and yy, is the matrix viscosity. 
Buckling is assumed to happen if the critical force, defined in Eq. (2.36) is 
reached, so 
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FE nl)" 
Ab = Ed (—2[Dijla:9;) 


= Br" (—2[D,;]]a:9;) 
>1, 


(2.38) 


with BP" being the dimensionless buckling index and DA is the normed 
strain rate defined as ID; jl = D;;/V. The index n corresponds to the possible 
fiber length Lf,. According to Figure 2.12 itis -1 < -2]|D;;]q:q; < 1 for all 
possible orientations in case of pure shear. Hence there is no state satisfying 
Eq. (2.38) if BP" < 1. Furthermore, Phelps et al. assume -2| D; qiq; =1in 
regions where the fibers are under compression. Regarding Figure 2.12, this 
means that the areas with x3 > 0 and x, < 0 or x3 < 0 and x, > 0 are fully 
and homogenous blue. Within this assumption, the remaining buckling 
criterium is BP" > 1 and no fiber orientations must be determined. 


Based on numerical experiments Phelps et al. [70] define a breakage proba- 
bility PP" for every possible fiber length as 


T 0 for BDU < 1 
P= Voy (1 - exp(1- Bb") for BDU > 1’ we 
with the breakage coefficient Cpr. Within [70] it is clearly mentioned, that 
(1 — exp(1 — Byu)) is a fitting approach, and better functions may be found, 
although the numerical experiments show a low sensitivity towards the exact 
formulation of this function. 


In a constitutive breakage model, fibers must create new, shorter fibers when 
breaking. Therefore, Phelps et al. [70] introduce a child generation rate 
defined as 


Lf 
RPE, = bt PDF (Sr) (2.40) 


being a Gaussian probability density function (PDF) with standard deviation 
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Spr and mean value 1/2, so fibers break most likely in the middle. The 
scaling factor @>F is defined to fulfill 


> Rim = 2P", (2.41) 


n 


since every break creates two children fibers. Furthermore, it is 
Ram = Rom-nyms (2.42) 


because the length of one child and the parent fiber defines the length of the 
other child fiber. For example, if one child of a 1.0 mm fiber is 0.7 mm, the 
other one must be 0.3 mm and hence have the same probability. The last 
restriction to the child generation rate is 


RE, =0 forn>m, (2.43) 


otherwise recombining fibers to create longer ones would be possible, which 
is unphysical. 


Finally, the evolution equation of fiber lengths is given by 


ON; ant SE : 
tu mem +) Rt Ni, (2.44) 


where Nf is the number of fibers with length Lf, within the control volume. 


In summary, the approach of Phelps et al. [70] is able to capture the evolution 
of an arbitrary number of different fiber lengths within an injection molding 
simulation, considering process conditions. Furthermore, only a small 
amount of material information is needed. Physical parameters such as, 
matrix viscosity, elastic modulus and diameter of the fiber are usually known. 
In addition, the three empirical parameters Ç, Cpr and Spr must be deter- 
mined. Although this approach is well suited to determine the fiber length 
distribution with respect to process parameters, there are disadvantages. One 
is the simplification —2[D;;lla:9; = 1 and another one is the need of an 
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empirical parameter for force calculation. These two aspects will be ad- 
dressed with novel approaches in Section 3.4 within this work. 


2.3 Simulation of process steps following 
mold filling 


The simulation of following process steps is indispensable to predict shrink- 
age, residual stress and hence the final part geometry. The range of complexi- 
ty and numerical effort of such simulations form a wide spectrum. The 
material model can be isotropic, linear elastic up to anisotropic, time- 
dependent and chemo-thermal viscoelastic. The simulation model itself may 
be only the filled cavity or the complete tool with integrated cooling/heating 
system. Several research works investigate in modeling with linear elastic 
[73], viscoelastic [74] and thermo-viscoelastic approaches [17,75-78]. 


In general, the material flow during the holding pressure phase is neglected 
and the simulations are pure solid mechanics. Besides the stress-strain equi- 
librium, the thermal equation (Eq. (2.3)) still needs to be solved, although the 
viscous heating can be neglected due to the low strain rates, compared to the 
mold filling phase. Furthermore, evolution equations for chemical reactions 
and phase change are still considered, since the majority of such processes 
take part within these process steps and they have significant influence on the 
material’s mechanical and thermal behavior. 


2.3.1 Matrix material modeling 


The mechanical behavior of the matrix is an overlay of material aspects such 
as the exothermal curing reaction, glass transission temperature, material 
history, but also molecular aspects such as polymer-chain geometry. Howev- 
er, the latter are not considered in macroscopic simulations. Besides material 
behavior, also the temperature field is significant, since it co-defines the 
material temperature distribution. Therefore, a wide range of models with 
different complexity and attendant experimental effort for parameter deter- 
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mination exist. A good overview of polymer viscoelasticity is given by 
Bergström [79] and in the special case of epoxies with fiber reinforcement by 
Bernath [14]. 


2.3.1.1 Viscoelastic modeling 


A viscoelastic material model is able to describe the material behavior within 
a wide range of frequencies. In addition, polymer phenomena such as relax- 
ing or creep can be modeled. In a general form the stress tensor is given by 


dex() 
dé 


t 
Oi = | Gje(lt - $) dé, (2.45) 
0 
with strain tensor €,;, time integration variable € and relaxation modulus 
tensor G;jr;, which is also often represented by a scalar, since the matrix 
behaves isotropic. Descriptions for G;;,; in case of injection molding warpage 
modeling can be found in [75,77] or by Choi and Im [17]. Unfortunately, all 
these studies focus on thermoplastic materials. According to the author’s 
knowledge during writing this thesis, no work addressing warpage simulation 
for reactive injection molding with visco-elastic phenolics has been published 
so far. Although, Mirabedini et al. [80] investigate in phenolic/rubber com- 
pounds, the viscoelastic behavior is dominated by the rubber part, and the 
phenolic is already completely cured within the studies. Nevertheless, viscoe- 
lastic modeling of continuous fiber reinforced epoxies is investigated in 
several studies and the general modeling aspects are transferable, since 
epoxies show a more comparable behavior to phenolics, than thermoplastics 
to phenolics. For thermoset materials, the relaxation modulus is often approx- 
imated with a Prony series [81] 


N 
t 
Gig (tT, c) = Gir (c) + Z Gijkip (C) exp (sma) (2.46) 
p= 


with Gj; being the elastic modulus in the rubbery state and ©, as actual 
relaxation time, depending on temperature and state of cure. The suitability of 


Prony series to model thermoset viscoelasticity is shown in [82-84]. Since 
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the relaxation time strongly depends on the material state [14,83], they are 
often approximated with a shift function such as 


8,(T,c) = a,(T,c) + OF (Tref, Cref), (2.47) 


where ort are relaxation times referring to a reference state based on a 
master curve determined with experimental measurements. The shift-factor 
a, can be modeled WLF based or with an Arrhenius-type function [14]. 


2.3.1.2 CHILE model 


Although viscoelastic approaches represent the material adequately, a high 
experimental effort is necessary to characterize the material and feed the 
model. Therefore, simplified approaches with accompanying restrictions 
exist. One example is the class of the cure hardening instantaneously linear 
elastic (CHILE) models. In the line of simplifications between fully viscoe- 
lastic and CHILE, also path-dependent models exist, but since they are not 
relevant within this work they are not further discussed. More information 
about path-dependent models is given by Svenberg and Holmberg [85] as 
well as Bernath [14]. 


Compared to viscoelastic approaches, the material history is neglected by 
CHILE models, so 


i dek (ë) 
aj = | Eğ (T, c) T dé, (2.48) 
0 
with elastic stiffness tensor defined as 
Eiko AT, 
Te(c) -T — Ta 
Egal) = 4 Efa + — (Ega Ega) Mr. (2.49) 
c2 ci 
Fey AT; 


Here it is AT, € T,(c) — T < Ta, ATs € Tg(c)-T > Tez and AT, < AT, < 


. .m C 
AT}. To, and Tez are material specific parameters. Efiki E; pa and Eiri are the 
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instantaneous elastic moduli in pure rubbery, pure glassy and actual material 
state. 


One disadvantage of the CHILE approach is the non-release of frozen-in 
stress, after reheating above Tg + T.ı, although there is a change of material 
stiffness. For normal process conditions, this aspect is uncritical, since the 
material is heated up to mold temperature once and cools down to room 
temperature afterwards. Both process steps include monotonous temperature 
changes. During heating, the degree of cure has not reached the point of 
gelation in most cases. Hence, no residual stress can be frozen in. In the 
cooling step, the change from rubbery to glassy state should only happen 
once, so the frozen stresses cannot release. 


The approach can be further improved by assuming Ejjx, and Ei not to be 
constant. Bogetti and Gillepie [86] describe Eiri by an interpolation between 
the values of fully cured and fully uncured, as function of the degree of cure. 


Adolf and Martin [87] describe the evolution based on the power law. But, 
due to a lack of experimental data, Efx, will be assumed to be constant 


within this study. According to Bernath [14] Ek 
constant, since the state of cure has neglectable effect on small-strain proper- 
ties, being relevant for modeling residual stresses and warpage. Furthermore, 
Bernath verifies the CHILE model to be well suited to model the residual 


stresses und warpage in case of transverse isotropic material behavior with 


ı is also assumed to be 


epoxy matrix in the RTM process. 


2.3.2 Homogenization and mapping 


2.3.2.1 General procedure 


As mentioned, the parameters have crucial impact on the material during 
processing, and therefore on the final properties of the part. Hence it is 
important to describe the material adequately and transform corresponding 
attributes from the process to the structural simulation, which are often not 
performed within the same software. 
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In case of structural simulation of FRPs, an additional material homogeniza- 
tion, based on fiber volume fraction and fiber orientation distribution, is often 
performed, because, similar to the process simulation, the numerical effort to 
simulate individual fibers is too high. Therefore, the homogenized FRP is 
described by an isotropic, transverse isotropic or orthotropic material model. 
The principle of homogenization is shown in Figure 2.13a. 


Since the meshes at start and end of the data transformation (also known as 
source and target mesh) are not identical, values must be interpolated. This 
procedure is called mapping, applied for example on temperature and curing 
field, but also fiber orientation tensors, relevant for the mechanical behavior. 


(a) (b) 


Figure 2.13: Schematic illustration of homogenization (a). Highlighting the influence of different 
meshes for change of values during mapping (b). Colors indicate different 


values of attributes 


Within the mapping process information may vanish or blur, as schematically 
shown in Figure 2.13b. The differences of source and target mesh have a 
significant influence within the mapping, as indicated in the intermediate 
state in Figure 2.13b, although this influence may become irrelevant, as 
indicated on the right side. While the mapping of scalars and vectors is 
explicit, this is not the case for tensors. 
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2.3.2.2 Homogenization schemes for FRPs 


The aim of homogenization is to calculate effective material parameters with 
accompanying minimal effort for experimental data or empirical model 
parameters. Especially in case of FRPs, homogenized material models are 
often used, since on the one hand, the numerical effort to simulate single 
fibers and matrix separately is too high and on the other hand the experi- 
mental effort to get all necessary material parameters for every material 
component is also unrealistic. 


The literature offers several homogenization schemes for different material 
groups and load scenarios. In general, homogenization schemes are separated 
in mean field and full field approaches. The latter are based on representative 
or statistical micro models to determine the mechanical answer of a macro- 
scopic load, as for example presented by Müller et al. [88,89]. Due to the 
experimental and numerical effort, these approaches will not be further 
discussed. Mean field methods can be divided in bounding and estimating 
methods. The most prominent boundary methods may be the upper bound by 
Voigt [90] in combination with the lower bound by Reuss [91]. Since these 
approaches only take fiber volume content into account, but no aspect ratio or 
fiber orientation distribution, the usage for discontinuous reinforced material 
is not suitable. Another often used approach for second order bounds is 
provided by Hashin and Shtrikman [92-94]. 


Additionally, some often used, well suited and more detailed approaches are 
the approach by Mori and Tanaka [95] and by Tandon and Weng [96]. The 
approach of Tandon and Weng considers fiber volume content and aspect 
ratio and does not need further empirical parameters. Based on the assump- 
tion of aligned spheroidal inclusions, a transverse isotropic stiffness tensor is 
determined for equal strain. 


Thermal properties can also be homogenized with different approaches. The 
specific heat can be simply volume averaged, since it is a scalar and isotropic 
property. The thermal conductivity can be determined with the approach 
given by Clayton [32], representing an extension of the Maxwell approach 
for dilute suspension and valid for higher concentrations. For the coefficient 
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of thermal expansion, Schapery presents an energy based approach [97]. 
Another often used model is presented by Rosen and Hashin [98], determin- 
ing the effective properties by superposition of an isothermal problem with 
surface displacement and an uniform temperature problem with homogenous 
boundary displacement. 


Nevertheless, these approaches only describe transverse isotropic behavior, 
which is insufficient in most cases for injection molded FRPs. The material 
properties must be further orientation averaged. Besides the orientation 
tensors, Advani and Tucker [51] present an orientation averaging scheme for 
tensor properties, based on a transverse isotropic microstructure. According 
to [51] a transverse isotropic second order Tensor T;; is fully described by 


Ti; = Bıqiq; + B2dij, (2.50) 


with the material specific constants B, and B,. Similar to the orientation 
tensor, described in Section 2.2.3.2 the orientation average of Tj; is given by 


(Tj) = | tuaova@oaa 


= B,(qiq;) + B2 (dij) 
= By Ajj + B26i;. 


(2.51) 


Here, (:) indicates the orientation averaging. Hence, the orientation average 
of a transverse isotropic tensor can be described by the orientation tensor and 
the corresponding unity tensor. 


A fourth order, transverse isotropic tensor is completely defined by five 
independent constants B,_s, as for example a unidirectional continuous FRP, 
where the mechanical behavior is represented by the five engineering con- 
stants. Therefore, the orientation averaging is given by 


(Tijen) = By Ajj 
+Bz(A;jöxı + Arıdij) 
+B; (Aikôjı + Andie + Adin + Aj Sir) 
+By(5ij5i) + Bs (Sixedj + Susi). 


(2.52) 
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The combination of the Tandon-Weng homogenization and the Advani- 
Tucker orientation averaging is successfully applied by Hohberg [99] for 
SMC material and also used within this work. 


2.3.2.3 Mapping schemes 


Several commercially available mapping tools exist to transform data be- 
tween a process simulation and a structural analysis. In case of fiber rein- 
forced polymers, these are for example MPCCI MapLib [100,101] or Digi- 
mat [102]. Further approaches are presented in literature, for example by 
Reclusado and Nagasawa [103], mapping fiber orientation based on only the 
first eigenvalue with corresponding eigenvector. MPCCI MapLib has been 
successfully used for SMC material by Hohberg [99] and Görthofer et al. 
[104]. The principal of these mapping processes is always the same and 
shown in Figure 2.14. 


Geometry and RS 
attributes of 


source mesh [WE 
_ OO r 
algorithm 


Figure 2.14: Schematic illustration of mapping process. Colors indicate different values of 
attributes. Due to source mesh, every element in the target mesh shows a dif- 
ferent color / value 


MPCCI MapLib maps independent of the software, which created the source 
and target mesh. Therefore, the mapping is purely based on geometry and the 
results and meshes. In case of mapping between different software packages, 
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the meshes must be converted in a neutral data structure [100], being for 
example VTK within this work. In case of tensors, MapLib interpolates each 
tensor component independently. This method may lead to a loss of shape of 
the tensor. Nevertheless, the only tensor mapped within this work is the 
orientation tensor, being positive definite and having a normed diagonal. 
Therefore, this mapping strategy is acceptable. 
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3.1 Considered interactions in process and 
material 


A comprehensive process simulation of FRP injection molding includes more 
than solving the Navier-Stokes equation with standard material models, as 
described in the Sections 2.2.1.1, 2.2.2 and 2.2.3. Figure 3.1 shows all con- 
sidered interactions within this work. At first, interactions between the FRP 
and the tool must be considered. This aspect is solved via boundary condi- 
tions, for example the tool temperature as boundary condition for FRP’s 
temperature at the boundary faces. Such interactions are already considered 
in a state-of-the-art single-phase and isotropic injection molding simulation. 


Due to the single-phase approach, the air in the mold is neglected in state-of- 
the-art simulations. Within this work, the air is considered as separate phase 
in the mold, leading to a multi-phase approach with two immiscible fluids for 
the simulations. This leads to two more interaction pairs, which need to be 
mentioned. One is, similar to the FRP, the interactions between the air and 
the tool, which are also solved via the boundary conditions. The other one is 
the interaction between FRP and air, which is solved within the multi-phase 
approach as described in Section 3.2 and [25,28]. 


Furthermore, internal material interactions are regarded within this work. The 
interactions between fibers and matrix are separated, so the influence from 
fibers on the matrix and vice versa are modeled individually and with differ- 
ent approaches. 


Due to their non-spherical shape, the fibers induce stresses in the matrix in an 
anisotropic way. Therefore, the FRP’s viscosity is modeled by a fourth order 
viscosity tensor, depending on non-Newtonian matrix viscosity, fiber orienta- 
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tion, length and volume fraction. The viscosity tensor and the underlying 
micro-mechanical models are described in Section and [26,38,105]. 


Tool and FRP 
> Boundary Conditions 


FRP and air 
> Multiphase approach 
(Section 3.2) 


Tool and air 
> Boundary Conditions 


Yj 
Z 


Fibers on matrix 
> Viscosity tensor 
(Section 3.3) 


Matrix on fibers 
>Hydrodynamic forces 
(Section 3.4.1) 


Fibers on fibers 
>Fiber-fiber forces 
(Section 3.4.2) 


Figure 3.1: Macroscopic interactions between FRP (green), air (light blue) and tool (shaded) with 
corresponding simulation aspects. Material internal interactions between fibers 
(blue) and matrix (yellow) with corresponding simulation approaches 


On the opposite, the relative velocity between matrix and fibers induces 
hydrodynamic forces in the fibers, leading to fiber re-orientation and may 
also lead to fiber breakage. The hydrodynamic forces are calculated to 
support fiber breakage simulations, depending on matrix viscosity, flow field 
and fibers. The corresponding methods are described in Section 3.4.1 and 
[37,106]. 


Besides the matrix, fibers interact with other fibers. These interactions induce 
lubrication, friction and normal forces in the fibers. Hence, these forces are 
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approximated as function of matrix viscosity, fiber orientation, geometry and 
volume fraction, as described in Section 3.4.2 and [106-110]. 


To keep the numerical effort in an acceptable range to perform a macroscopic 
process simulation of a whole part, the fiber orientation is described with 
state-of-the-art models and orientation tensors, described in Section 2.2.3. 
Hence, for the hydrodynamic forces on fibers, only the portion in fiber 
direction is considered, leading to buckling and breaking. The non-fiber- 
direction component of the forces leads to re-orientation and would therefore 
be an overshoot, because the orientation is calculated separately. 


3.2 Multiphase approach! 


3.2.1 Volume-of-fluid method 


A multi-phase approach requires a methodology to differentiate between the 
phases and separate material properties in the calculations. Many different 
approaches to separate two, or more, phases in a CFD simulation have been 
published and can be found in several fluid dynamic books. Within this work, 
the FRP and air are considered as homogeneous fiber matrix suspension and 
gas mixture. Therefore, the multi-phase approach distinguishes only between 
two phases, which are immiscible. A well-suited approach for these condi- 
tions is the volume-of-fluid method (VOF) with one VOF-factor a € [0,1]. 


Within this work, a = 0 equals pure air and a = 1 pure FRP. Cells with a € 
(0,1) are partially filled. Effective material specific properties can be calcula- 
ted by volume averaging 


"eff a|: [FRP + (1 u @)| lair: (3.1) 


! The methods presented in this Section are published in [25]. 
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3.2.2 Phase-dependent boundary condition 


Due to the presence of air within the mold, a phase-dependent boundary 
condition must be formulated, the way that air can leave the mold but FRP 
not. The boundary condition is defined as function of @ to consider the phase- 
dependence and applies for U;. The air is assumed to leave the mold without 
a change of state, so 


au, (2 9 0 
ae 0 0 O} fora =O, (3.2) 
% \o 0 0 


j 
being a Neumann boundary condition named zeroGradient. Furthermore, the 
boundary condition should act like a wall for the FRP resulting in a Dirichlet 
boundary condition with 


U; = (000) fora = 1, (3.3) 


named U,ero. Hence a no-slip boundary condition is applied, which is a 
simplification for RIM, but creates good results with less need of model 
parameters (see Section 2.2.1.5). 


There are different options to define the condition for a € (0,1). One is to 
define a fixed switchover for a = 1, leading to a loss of material caused by 
outflow, while cells are partially filled. To reduce this loss, the switchover 
can be defined at any point œ < 1, which leads to in-mold air, unable to 
leave. Within this work, a compromise between these two cases is chosen. 
Similar to material properties and according to [25,26] the two boundary 
conditions zeroGradient and Uzero are interpolated with a. The final velocity 
boundary condition is described schematically with 


Ubound = @U zero + (1 — a)zeroGradient. (3.4) 


Numerical studies and a verification of the boundary condition are presented 
in Section 4.1.1. To minimize the material loss caused by this formulation, 
the condition is only applied for a specific region of the mold as shown in 
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Figure 3.2. This region, called outlet, is defined at the end of the flow path, 
where the venting slots are positioned in the real process. The other boundary 
areas (besides the inlet) are impermeable walls, so neither FRP nor air may 
leave the mold. 


e inlet 
wall 
e outlet 


Figure 3.2: Different boundary regions for an arbitrary mold filling model. Inlet in blue, imper- 
meable walls in grey and outlet with phase-dependent boundary condition in 
red 


3.3 Anisotropic viscosity modeling? 


3.3.1 Theory 


3.3.1.1 Transversely isotropic fluids 


The real viscosity of a FRP is a complex combination of temperature, chemi- 
cal reaction, shearing and, due to presence of fibers, anisotropy. Within this 
Section the focus is on description of anisotropic viscosity due to fibers. 


Since fibers are non-spherical, stress builds up along the surface in an aniso- 
tropic way. This effect is similar to FRPs in solid mechanics and hence, the 
simplest anisotropic behavior is transversely isotropic. This behavior is 
present if all fibers are aligned parallel in one direction (cf. Figure 2.11b and 
Figure 3.3). Gibson [43] describes a transversely fluid analogy to solid 
mechanics by 


? The methods presented in this Section are published in [26]. 
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&j = Si Cid 
Syn $1122 $1122 0 0 0 
S2222 $2233 0 0 0 (3.5) 
transverse = S2222 0 0 0 
gs 453322 — S1111 0 0 
sym. S1212 0 
S1212 
with S;jxı as forth order fluidity tensor. 
Additionally, Gibson assumes incompressibility, so it is 
£1 + £22 + €33 = 0, (3.6) 
leading to 
(022 + 033)(S2222 + S1122 + S2233) 3T) 


= —011 ($1111 + 254122). 


Since Eq. (3.7) must be valid for any arbitrary stress state, two more condi- 
tions are given by 


(S2222 + $1122 + S2233) = 0, 
and (3.8) 
(S1111 + 251122) =0, 


hence the number of independent material constants reduces to three and the 
fluidity tensor can be rewritten as 


çtransverseincomp en 
ijkl = 


51111 51111 
Sı111tS2323 Sı111-52323 0 0 0 
4 4 
Suuts3 o o 0 (3.9) 
4 ’ 
$2323 0 0 
sym. Sp» 0 
S12 12 
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Within [43], these three constants are independent of fiber attributes. 


3.3.1.2 Dependence on fiber attributes of transversely isotropic 
fluids 


A solution to describe the three constants of a transversely fluid as function 
of fiber attributes is given by Pipes et al. [105]. Besides the full alignment 
and the incompressibility Pipes et al. define a few more assumptions. The 
fiber ends are touching and the fibers are arranged in square or hexagonal 
packaging. Fiber ends in one row are next to fiber centers in neighbor rows, 
also in hexagonal packaging, which is a simplification. All fibers have the 
same diameter df , length Lẹ and neighbor distance h;. The rows have an 
offset of Le/2. The variation of velocity in fiber direction is linear (U, = 

1X1 + const.), so €;, = const. Fibers move with the velocity of their 
center. A schematic visualization of such an assembly with resulting shearing 
is given in Figure 3.3. 


Uy 2 
dr Lr/2 
D 
— —_ —_ =N h 
x = = I 
| 2 D 


Figure 3.3: Aligned fibers in a fluid. Shearing due to U12 > U11. 


Based on the assumptions and information in Figure 3.3 it is 
U12 — U11 = &ıLe/2 (3.10) 
and 


12 =Ý = &ı1Lı/(2(he - dp). (3.11) 
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For a fiber formation as shown in Figure 3.3, the fiber volume fraction is 
given by 


d 2 
Di = Dmax (=) (3.12) 


with Ọ nax being the maximum possible fiber volume fraction, which is 
n/N12 for hexagonal and 7/4 for square packing. Since hexagonal packing 
represents a stable equilibrium it is the more realistic case. Therefore, hexag- 
onal packaging is assumed within this work. 


According to [105] the average force F acting in fiber direction at the fiber 
midpoint is given by 


a 20 (nd? 
F= —— |. 3.13 
a ( 7 (3.13) 


On the opposite, the force due to shear stress T,2 on half of the fiber surface 
is 


A tT Led t Led 
e Te aw 


Building up the force equilibrium leads to 


p “Oji. *) O11 
= —|— | = — . wl 
MY Dr (i Pere G 5) 


Finally, the combination of Eq. (3.11), (3.12) and (3.15) yields to 


un O11 a sii V P/P max ) 2 (3.16) 
See f: : 
E11 2 1- Vy P/P mnax 


The axial elongational viscosity 71, is identical to the entry 71111 in the 
fourth order viscosity tensor for a complete micro-mechanical approach. 


711 


Nevertheless, later in this work the formulation of 71111 will be defined 
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differently due to homogenization based on this micro mechanical approach 
(see Section 3.3.2). Therefore, the notation n41 is chosen. 


Similar to 7, the axial shear viscosity 712 and transverse shear viscosity 723 
can be determined, as shown in [105]. The results are 


u NM 2- D/P max 
N12 = Der Tora: (3.17) 
L= y P/P nax 
and 
n 
123 = - (3.18) 


(Mi u Dax)? 


With these three parameters (N11, N12 and 723) it is possible to fully describe 
the viscous behavior of a transversely fluid with fiber arrangement as shown 
in Figure 3.3. The matrix viscosity can be assumed constant, as in [38], or be 
calculated with an viscosity model such as in [26]. Within this work, the 
matrix viscosity is calculated with the Castro-Macosko model described in 
Section 2.2.2.3, if not explicitly mentioned differently. 


Pipes et al. published further investigations to describe the material behavior 
on microscopic scale for a fiber-fiber overlap + L¢/2 [111] or within a 
power-law matrix [112]. Nevertheless, these further approaches are not 
suitable for a macroscopic process simulation with orientation tensors, since 
too much information, for example the overlap length of individual fibers, is 
not defined. 


3.3.2 Viscosity tensor of discontinuous FRPs 


3.3.2.1 Viscosity tensor of transversely isotropic FRPs 


Based on the assumptions and formulations of Section 3.3.1 and by setting 
S1111 = 1/111, $2323 = 1/23 and S1212 = 1/Nı>2 the viscosity would be 
defined as inverse of the fluidity tensor (Eq. (3.9)), but, due to the incom- 
pressibility assumption, the fluidity tensor is singular. Sommer et al. [38] 
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present a solution to formulate a fourth order viscosity tensor for the assump- 
tion of incompressibility. In a first step, a pseudoinverse viscosity tensor for a 
transversely isotropic fluid is built. The methodology is given by Loredo and 
Klöcker [113], determining the stiffness tensor of an incompressible solid 
material. The pseudoinverse is defined as 


-1 


1 1 
Sji = (Sina + (3) 65x) = (3) ijn: (3.19) 


Applying this for the fluidity tensor as defined in Eq. (3.9) and further replac- 
ing $1111 = 1/M11, S2323 = 1/23 and S1212 = 1/12 leads to 


Amıı 211 — 2711 0 0 0 

111 + M2311 — 9N23 O 0 0 

transverse _ 1 N11 + 923 0 0 0 
Nijkl 5 nog 0 0 (3.20) 

sym. 12 0 

12 


3.3.2.2 Orientation averaging of the viscosity tensor 


The tensor given by Eq. (3.20) is only valid for transversely isotropic, mean- 
ing fully aligned, materials. However, this state is never reached in a real 
process, although it might be a good approximation near the surface and in 
thin part regions. To close this gap and apply the tensor also for less orientat- 
ed regions, the orientation averaging by Advani and Tucker, introduced in 
Section 2.3.2.2 (Eq. (2.52)) is used, similar to [26,38]. Finally, the orientation 
averaged fourth order viscosity tensor is defined as 


(nije) = (mı — 412 + N23) Aijx 
n 
ag (- a + n23) (Aijôki + Aridi) 
Aix Sj, + 
+ = 
(12 = Nea) ag Anda Ade 

M11 
+ (> = n23) ($;;6x1) 
+23 (Sirßjı T ud). 


(3.21) 
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This formulation fulfills Eq. (3.9) and generates similar results compared to 
the formulations of Ericksen [114] and Hinch and Leal [115]. The viscosity is 
calculated anisotropic, depending on fiber orientation, fiber volume fraction, 
fiber aspect ratio and matrix viscosity, represented by the three parameters 
711, N12 and 23 given in Eq. (3.16)-(3.18). 


3.3.2.3 Extension to multi-phase modeling 


The work of Sommer et al. [38] focusses on modeling a single-phase com- 
pression molding process. The approach must be extended to be applicable 
for a multi-phase simulation with FRP and air, as described in [26]. The 
viscosity is a material specific variable, interpolated with the VOF method 


(cf. Section 3.2.1). Hence the effective viscosity ne is 


nif = a nih + A- a) ni (3.22) 
with 
4/3 -2/3 -2/3 0 0 0 
4/3 -2/3 0 0 0 
i 4/3 (0) 0 0 
nait, = Mar / E a er 
sym. 1 0 
1 


where Nair is the scalar viscosity of air, assumed to be constant within this 
work. The fourth order viscosity tensor of air is isotropic and the formulation 
fulfills the condition 


nei Di = 
1 dU, (3.24) 
Mair (Di; ETA u) 


so the results for the viscous stress tensor are identical for the scalar and 
fourth order tensor formulation (see Section 2.2.1.1, Eq. (2.4) and Eq. (2.5)). 
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3.3.2.4 Implementation in OpenFOAM 


The momentum equation (Eq. (2.2)) is solved implicitly within OpenFOAM 
4.1. Since the viscous stress tensor directly depends on the velocity gradient, 
Eq. (2.4) is used during solving, so the viscous stress can change and the 
equilibrium is reached. This term cannot be directly replaced by nfk Der 
since the implicit solving algorithm is only defined for scalars, vectors and 
second order tensors in Cartesian coordinate systems. Nevertheless, implicit 
solving is more stable and robust. To gain more stability during solving, the 
effective viscosity tensor is split into an isotropic and anisotropic part so that 


eff _ „effiso eff,aniso 
Niji = Nike FN > (3.25) 


as presented in [26]. The isotropic part can be represented by a scalar value 
and hence be solved implicitly without loss of information. According to 
Bertöti and Bohlke [116] the effective scalar shear viscosity of an anisotropic 
fluid is 


1 
Neffiso = To Mile Pör (3.26) 


where Pd is the second projector tensor of fourth order. The corresponding 


fourth order tensor is defined similar to the viscosity of air as 


4/3 -2/3 -2/33 0 0 0 
4/3 23 0 0 0 
i 4/3 0 0 0 
i Sas hr 2) 
sym. 1 0 
1 
The anisotropic part can be determined by 7 aa = ne — ne °°, Finally, 


the viscous stress is calculated with 
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„eff; = eff,iso eff,aniso 
Tij = NijkıDrı = (niiki + Nijki Dx 


1 dU, 
ff,aniso 
= Mer; (Di -55 ôy) + E ASS e (3.28) 
Neffiso ij 3 Ox, ij Nijki kl 
solved implicitly solved explicitly 


3.4 Fiber forces and breakage modeling? 


3.4.1 Hydrodynamic forces from fluid on fibers 


The anisotropic viscosity tensor, presented in Section 3.3 represents the fiber 
to fluid part in the stress equilibrium, but not the fluid on fiber part. The 
matrix introduces forces in the fibers, due to relative velocity, leading to re- 
orientation or in some cases to fiber damage. This fluid dynamic phenome- 
non is known as hydrodynamic forces. In principal the effect is similar to 
aerodynamics and therefore, the same approaches can be adapted. 


Since fibers are non-spherical bodies, the hydrodynamic force Ee is a 


combination of drag force FA and lift force FË. The first appears whenever a 
fluid has a relative velocity and contact to a rigid body, while the latter occurs 
due the inhomogeneous stress distribution on the surface, caused by the non- 
spherical geometry [20]. 


The literature offers several studies on the modeling of hydrodynamic forces 
on fibers, such as for example work of Phan-Thien et al. [117], Lindström 
and Uesaka [118], Dinh and Amstrong [72] and Meyer et al. [37]. Within this 
work, the approach of Meyer et al. is adapted, since it has proven to be well 
suited for modeling forces on glass fibers, or fiber bundles, in a macroscopic 
process simulation, although Meyer et al. focus on compression molding 
simulations. Furthermore, the approach provides the possibility to distinguish 


3 The methods for fiber force modeling presented in this Section are pub- 
lished in [106]. 
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between drag and lift force, which is not directly possible in other approach- 
es. 


3.4.1.1 Modeling ofhydrodynamic drag and lift force 


The drag force of a sphere is given by Stokes law, so 
Ff = 3nnmd.pAU;, (3.29) 


with dsp being the diameter of the sphere and AU; the relative velocity. 
According to Batchelor [20], the absolute hydrodynamic resistance is propor- 
tional to Nudsp|lAU;|| and independent of the body’s shape, if the flow is 
Newtonian, incompressible and the internal forces are neglectable compared 
to the viscous forces so Re « 1. Hence, Meyer et al. define an equivalent 
diameter 


Asp = kal b) ds, (3.30) 


with the dimensionless correction coefficient kg, depending on fiber aspect 
ratio and relative angle between fiber and relative velocity &. 


Besides the drag force, the lift force should also be considered. Meyer et al. 
[37] present a similar approach using a dimensionless correction coefficient 
ky. Of course, the lift force acts not in direction of the relative velocity, but 
perpendicular to the relative velocity and the fiber axis. Hence, the lift force 
direction is defined as 


pi = (qi X TAU; ]) x [AU]. (3.31) 


Similar to the drag force, but acting in a different direction, the lift force is 
then defined by 


FË = 3nnmki (re &)dellAU,|I [pi]. (3.32) 


Finally, the complete hydrodynamic force on one fiber is given by 
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hyd i 
aa aes (3.33) 
= 3nmde(kaAU; + kylAU;lILpil). 


The [-] brackets represent a normed vector, defined as [x;] = x;/IIxill- 


3.4.1.2 Definition of correction coefficients ka and kj; 


Within [37], kg and kj; are approximated, based on numerical experiments 
computing the stress on single fibers in a simple shear flow with different 
orientations and aspect ratios. Since the hydrodynamic force is the result of 
an interaction between fluid and fiber on the fiber surface, it is also given by 


R = | onda (3.34) 


l 


with Aç describing the fiber surface and n; being the corresponding normal 
vector. Here, the fiber ends are neglected. Meyer et al. [37] performed nu- 
merical experiments with a fixed fiber within a flow field with known viscos- 
ity and velocity. Combining Eq. (3.33) and (3.34) leads to 


1 
A ET y) . 
a it ER en 
and 
k Er rn, | dA 3.36 
i= 3rd lau | Mer me) 


where the index i and j indicate the components vertical and horizontal on 
the surface here. Meyer et al. solve Eq. (3.35) and Eq. (3.36) for AU; = 
(1 0 O)mm/s and ¢ € [0°,90°] as well as different aspect ratios. The 
angle interval represents all possible angles between fiber and relative veloci- 
ty, since a fiber has no front or back end. Within [37] numerical fits are 
performed to approximate kg and kj; as function of aspect ratio and angle. 
The results are given by 


ka = 1 — anya (re — 1)cos(2$) + Bnyalrr — 1) G37) 
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and 
kii = Anya; — 1)sin(2¢), (3.38) 


with hya = 0.09 and Phya = 0.3125. By these approximations kg shows a 
steady increase for rp > 1 and for an increase of @, while it is constant for 
rr= 1 , applying [20] in case of spherical inclusions. Whereas kj; is zero for 
rp = 10rd = {0° 90°} and has a maximum for & = 45°, since fibers have 
a symmetric surface. 


3.4.1.3 Application of hydrodynamic force modeling in 
homogenized material 


The work of Meyer et al. [37] includes a material model approach with 
individual fibers (or fiber bundles) enabling the calculation of the hydrody- 
namic force as given by Eq. (3.33). Within a homogenized material three 
important aspects are unknown. The first one is that the relative velocity 
between fibers and fluid is unknown. Secondly, the angle between the fiber 
and AU; is unknown. Consequently and thirdly, the direction of the lift force 
pi is also unknown. In a first step, the relative velocity is approximated. 
Similar to [37] and [106], AU; can be approximated by 


W, 
AU; = w, (Uin — Uim), (3.39) 
n 


with Wyrm = exp (-Od3m)/ C) and Wn = YinWnm- The indices n and 
m indicate different cells and dpn is the distance between the centers of these 
cells. Contrary to [37] Uim is not the velocity of the bundle part within this 
work, but the velocity in cell m at the cell center and is assumed that all 
fibers within this cell have on average the same velocity as the cell center. 
This simplification is justified by homogenization and the fact that no fiber 
matrix separation is assumed [106]. Hence, AU; depends on the velocities 
Uin of the n considered neighbor cells. The value of n is chosen so 1/ 
NYinAnm = 1.5 Ls, but at least one generation of cell neighbors is always 
taken into account. 
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The next step is to approximate the angle between AU, and the fibers. With- 
out knowledge of single fiber orientation, the next meaningful approximation 
is to determine some kind of average angle. This angle can be defined in 
relation to a reference fiber. In structural analysis, the material is often 
described transverse isotropic with fibers along the first eigenvector, as 
described in Section 2.3.2 and [103]. Hence, the eigenvector can be interpret- 
ed as such a representative fiber. Since there are always three eigenvectors, 
there are three reference fibers with exact known orientation and the orienta- 
tion probability of the corresponding eigenvalue can be determined for every 
known second order orientation tensor. This assumption fits to the definition 
of the orientation tensor, since it is 


3 
Ay = | Wadaadq = X Anvinvin (8.40) 

n=1 
in case of three fibers. Here, Vin is the n-th eigenvector of A;j and An the 


corresponding eigenvalue. 


By knowing the orientation of the eigenvectors, the angle to the relative 
velocity is given by 


6 (752m) GB 41) 
= arccos B 
7 AU; || 
and the direction of the lift force by 
Pin = (Vin X [AU;]) x [AU;]. (3.42) 


Based on this information, Eq. (3.33) can be solved and me can be approx- 
imated in a discretized mesh with fiber orientation tensors [106]. Due to the 
dependence of @ and p; on the three eigenvectors, pe is calculated three 
times in every cell. Furthermore, kg and kj; depend on the fiber aspect ratio, 
so in case of fiber length distributions, al is determined three times for 
every possible fiber length. 
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3.4.2 Fiber-fiber interactions 


3.4.2.1 Approximation of fiber-fiber contact points 


During processing, fibers are in contact, which gets more significant for more 
highly filled suspensions. Whenever there is a contact between two fibers, 
contact forces act at these contacts. Since the forces appear on every contact 
point, in a first step, the number of these contact points needs to be deter- 
mined. 


Toll [119] describes that the average number of fiber centerlines, intersecting 
an average test fiber is exactly given by 


Na = neL2def + 1/4 mn¢Led?(g +1), (3.43) 


with ne being the number of fibers per volume and f and g are the scalar 
invariants of the fiber orientation distribution, given by 


f= [lat x aflvca caf aaa? (344) 
and 
g= | | lata? was) Wa) dasa? . (3.45) 


Here, the indices a and ß indicate different individual fibers. In a later study, 
Toll shows, that the average number of fiber contact points Npe per fiber is 
determined by replacing d; within Eq. (3.43) with 2d, [107], so 


Nec = 2nl2def + nneLed? (g + 1). (3.46) 


The fiber volume fraction &; is given by (1/4)nn;L;d?, so the contact points 
may also be calculated as function of the volume fraction by 


Nec = Em Oeref + 4O;(g + 1). (3.47) 
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Finally, the number of contact points is given, depending on fiber volume 
fraction, aspect ratio and orientation. Although this approach is a great 
advance in FRP process modeling, one disadvantage is the dependence of f 
and g on single fiber orientations, being an information, which is not availa- 
ble in a macroscopic process simulation. 


Therefore, Férec et al. [120] present an approach to approximate the number 
of contacts points, with information provided by the orientation tensors. 
Within the study, new rheological models for fiber orientation are developed, 


including a so-called interaction tensor bł, j» defined as 


= by = { | qrat la x qf lya Dya dqa.. (3.48) 


Due to the normed orientation vectors, it is 
f= b}. (3.49) 


Furthermore, Ch x qj Bl = =1- “es 


=| [Ge q% - agapatabatas 
waa} darda 
= | | ras — a asagarar a; ) 


w(qe)v (a) datda? 
= = Ai — ArjrıAkı- 


(3.50) 


Additionally, Férec et al. [120] introduce a fitting factor, so values for f 
calculated with Eq. (3.44) and Eq. (3.49) are identical in case of quasi- 
isotropic fiber orientation and the final interaction tensor is finally given by 


== ur (Aij z Ar Arı)- (3.51) 


According to Toll [107] the term 4;(g + 1) in Eq. (3.47) can be neglected 
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in case of high aspect ratios 77, only appearing in the first term, which be- 
comes dominant in this case. Férec et al. neglect this term too and no approx- 
imation for the scalar invariant g is given in [120]. Nevertheless, for high 
orientations f decreases, while g increases and becomes more important as it 
is also the case in short fiber materials. Within this work, the two novel 
approaches to determine f and g presented in [106] are used, also only 
depending on information, given by A;;. In this way, the average contact 
points of a fiber within a fiber network can also be determined adequately in 
case of high orientations and short fiber materials in a macroscopic simula- 
tion. 


In the first approach the eigenvectors Vin and eigenvalues A, of Aj; are used 


to determine f and g. Similar to Eq. (3.40) f and g can be reformulated as 


3 
pe > [vin X Yn|AnAm (3.52) 
n,m=1 
and 
3 
g= >. WinvmlAnAm. 3.53) 
n,m=1 


Since Ajj € sym. it is Vin L Vim, SO |VinVim| = 0 for n +m. By further 
regarding |v;| = 1, f and g simplify to 


f = "Jg (24142 + 24143 + 24243) (3.54) 

and 
g = àA AAgAg+A3A3. (3.55) 
As shown in [106] and 4.1.3.3 the calculation of f with the interaction tensor 


and Eq. (3.54) creates identical results. Therefore, the correction factor 37/8 
is also used for this approach, to fit the value in case of quasi-isotropic 
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orientation. In the case of full alignment, it is g = 1 for Eq. (3.45) and Eq. 
(3.55) and since this invariant is more important for more highly orientated 
states, no fitting factor is used. Eq. (3.54) and Eq. (3.55) show, that f and g 
can be approximated as function of the eigenvalues with corresponding 
factors in a polynomial way [106]. Hence, in a second approach, the invari- 
ants can also be approximated by 


f.9= >. MamÄnkm 3.56) 


with Mpm containing the polynomial coefficients and A, = (A, A 1), 
since the eigenvalues are normed and all information is provided by two 
eigenvalues. The entries of Mpm are determined in [106] and given in Table 
3.1. 


Table 3.1: Coefficients of polynomial fit of scalar invariants f and g [106] 


Entry of Mpm Approx. of f Approx. of g 
Mi, BD, 3.3011 
M32 -6.6744 0.4173 
M33 1.3475 1.5728 
Miz + M31 4.63897 4.2687 
Miz + M31 -4.5262 -3.8701 
M33 + M32 2.482 -1.9965 


3.4.2.2 Approximation of fiber contact forces 


The three forces, present at every fiber contact point, are the normal force F’, 
the friction force Fi" and the lubrication force pu The direction of the 
normal force depends on the orientations of the two fibers in contact. Since 
the reference fibers are represented by the eigenvectors within this work, the 
normal force cannot be approximated adequately due to A;, € sym. and vj, L 
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Vim for n # m. Nevertheless, the amount of the average normal force, as for 
example used in [121,108,109] and relevant for friction, can be given by 


32 
IF? || = zya Erde OF f’. (3.57) 


One benefit presented within [106] is the possibility to model the invariant g 
and not only f. Therefore, when modeling the normal force, and hence 
friction force, g can also be considered, as shown in [106]. According to Toll 
and Mänson [108] the average normal force within a volume is given by 


are | 
FP ll =-—_,, (3.58) 
: Nnodedr 
with the node force FR°@® and the number of nodes per volume n„ode defined 
as 
ddr 
Nnode = md? ay (3.59) 


where a; is the average node space. One node is defined by four fibers with 
three contact points, as schematically shown in Figure 3.4. By assuming the 
contact points to be homogenously distributed along the fibers and according 
to Figure 3.4 the average node space is given by 


_ 20 2L¢ „= 
dr = — = ; : 
‘Nic PeNice (3.60) 
with the specific number of contacts Nig (cf. Eq. (3.47)) defined as 
8 Nee 
Nico = /mrif +4g +1) = or (3.61) 


According to [108], FP°4® is given by 
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|jrpode || = = N 5 dor i (3.62) 


where Sp is the average nodal compliance, defined by 


— 243 


= ——. 3.63 
a MEd? ( ) 


Figure 3.4: Schematic illustration of one node as defined in [108], consisting of four fibers and 
with node space a. The contact points are assumed to be homogeneously dis- 


tributed along the fiber 


Solving Eq. (3.62) by considering Eq. (3.60) and Eq. (3.63) leads to 


dr 
node|l _ Erde pr “Nig 
a a 
0 
_ Edi bs Nio 
4014 


ddr 
(3.64) 


Further combining Eqs. (3.58), (3.59) and (3.64) the average normal force is 


ME de be Ngo 


FD = 
el =e 


(3.65) 


For Neg = 8% f/m, meaning neglecting the g-term of the contact points, 
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which is the assumption within [108], Eq. (3.65) and Eq. (3.57) are identical. 
Within this work and similar to [106] Eq. (3.65) is used for normal and hence 
friction force. 


Regarding the friction force, Coulomb friction is assumed in most cases and 
the force is given by 


FË = k„|FPII]JAUF], (3.66) 


with ker being the friction coefficient and Auf is the relative velocity be- 
tween the two fibers. Of course, the real relative velocity of the fibers is 
unknown due to the homogenization and the relative velocity AU; between 
matrix and fibers defined in Eq. (3.39), used for calculating the hydrodynam- 
ic force, is not valid in this case. The relative velocity of two fibers within 
one cell is given by 

AUF = Dijhf, (3.67) 
where hffis the distance vector between the fiber centers. Assuming the 
distance to be equal in every direction [106], Eq. (3.67) can be reformulated 


to 
nn hhf w 
Auf = BRS Dy = anf (3.68) 


J 


Since the friction force only depends on the normed relative velocity, nit does 
not need to be determined and the direction can be approximated by summing 
up the corresponding entries of the strain rate tensor. 


Regarding the lubrication force, the literature contains several approaches 
with different level of detail and complexity. Yamane et al. [122], Bounoua et 
al. [123] or Meyer et al. [110] present detailed approaches depending on 
single fiber orientations and therefore not suitable for a macroscopic ap- 
proach. Approaches suitable for a homogenized material are usually assumed 
to be linear proportional to a dimensionless coefficient kju. The linear de- 
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pendence is further on the relative velocity Auf and matrix viscosity, as for 
example presented by Servais et al. [109]. Djalili-Moghaddam and Toll [124] 
further consider the surface area, but the area is assumed to be constant and 
the dependence is represented by df. Within this work and identical to [106], 
the lubrication force is assumed to be proportional to the relative velocity, 
matrix viscosity, overlapping area and reciprocal proportional to the fiber 
distance, since it vanishes with rising distances. Finally, the lubrication force 
is given by 


k 
pe = nA 


Infnft /3 (3.69) 


= kun Ard0f, 
with Age being the average overlap area of two fibers within one cell. 


3.4.2.3 Approximation of fiber-fiber overlap area 


To approximate the fiber-fiber overlap area, it is assumed, that the fibers do 
not overlap near the ends [106]. Therefore, the overlap area is completely 
defined by the fiber diameter and the overlap angle w. Two different cases, 
depending on orientation state must be separated. The first one is w > Werit, 
where the overlap area is a parallelogram, as shown in Figure 3.5a and Figure 
3.5b. The second one is w < Wcrit, Where the area is a hexagon, as shown in 
Figure 3.5c. 


If itis > Writ the projected overlap area is given by 


Ag = an 
ft sin(w)’ 


(3.70) 
For realistic aspect ratios (rẹ > 10) and material orientation states, it is w > 
crit In Most processes. Nevertheless, the case w < w.rje can be reached, and 
must be considered. In this case, the assumption of small angles is an ade- 
quate approximation and the simplification sin(w) = w is valid, hence the 
critical angle is given by 
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w = Werit 


p 


\ 


A 
Z 


\ 


ZA 


(a) (b) (c) 


Figure 3.5: Two overlapping fibers with highlighted overlap area (green). Arbitrary angle w > 
Wcrit (a), critical angle w = werit (b) and over-critical angle w < werit (C). 
Red area is subtracted for calculation of overlap [106] 


Ort I => (3.71) 


The overlap area is given by subtracting the not overlapping areas from the 
complete projected surface. Regarding Figure 3.5c, the red areas are subtract- 
ed from the complete area and result is the green area. So, for W < Wecrit the 
overlap area is given by 


12 
Aff = deLe = J% (3.72) 


Combining Eq. (3.70) and (3.72) the projected overlap area is given by 


dê í " 
A or @ = Werit 
food SO) (3.73) 


12 
idiom is for sch, 
ff 4 crit 


The angle of the individual fibers is simply given by 
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w= cos”1(qq°). (3.74) 


Unfortunately, the eigenvectors are always perpendicular to each other, so 
computing the individual angle would not be meaningful or representative. 
Furthermore, the average w is unknown in the homogenized material, but as 
shown in [106] the averaged angle within one cell can be approximated with 
f and is given by 


w x —. (3.75) 


Hence, the averaged, projected overlap area of the fibers within one cell can 
be determined depending on fiber geometry and orientation state, with 
information provided by the second order orientation tensor, fiber length and 
diameter. 


3.4.3 Constitutive fiber breakage model 


3.4.3.1 Breaking criterium based on hydrodynamic forces 


As explained in Section 2.2.4.2, buckling is one of the major mechanisms for 
fiber breakage. In a first step, it is determined, whether the fibers within one 
cell are under pressure and buckling in general is possible. This state is 
reached if D;;q;q; < 0 (Eq. (2.34)). Similar to calculation of hydrodynamic 
forces, explained in Section 3.4.1, the eigenvectors of the orientation tensor 
are used as reference fibers, and the Jeffery orbit is reformulated to 


DijVinYjn < 0, (3.76) 


where, again the index n indicates the number of the eigenvector and n € 
{1,2,3}. A summation of the corresponding eigenvalues approximates the 
amount of fibers, that are able to buckle, within one cell. 


In a next step, it is checked, whether the force, acting on the fiber is high 
enough to cause buckling. This concerns only the eigenvectors fulfilling Eq. 
(3.76). Therefore, the hydrodynamic forces are used, as described in Section 
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3.4.1. For buckling, only the amount of force in fiber direction is relevant, 
which is given by the dot product of the force and the corresponding fiber 
direction. So, the relevant force is given by 


pt aq ah (3.77) 


and the buckling criterium is reached if 


(3.78) 


with E>" being the critical force for fibers with length Lf, and defined in Eq. 
(2.36). Similar to [69] and [70] transverse forces, which favor buckling, are 
not considered, which is a huge simplification due to a lack of information in 
an macroscopic and homogenized approach. 


3.4.3.2 Fiber length evolution 


The evolution of fiber lengths is strongly based on the model of Phelps et al. 
[70] explained in Section 2.2.4.3. The breaking probability is defined as 


a 0, caseA 
) a 
Pam m ds (1 = exp(1 a Bis i caseB’ (3.79) 


for every eigenvector Vin and fiber length Zf. The cases are 


hyd 


caseA = By <1 or DijVinYjn > 0 (3.80) 
and 
caseB = prya >1 and DijVinVjn < 0. (3.81) 


The final breakage probability for a specific length is weighted with the 
eigenvalues and given by 


Pra, Ber (3.82) 
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The subsequent fiber breakage simulation is similar to Phelps et al. [70]. The 
child generation rate is defined similar to Eq. (2.40) with the same re- 
strictions as described in Section 2.2.4.3 and [70]. The time evolution is 
given by Eq. (2.44), simulating the fiber length distribution by number of 
different possible fiber lengths in every cell. 


3.4.3.3 Model restrictions and separation to state-of-the-art 
modeling 


Within Section 3.4.3.1 only the hydrodynamic forces from Section 3.4.1 are 
mentioned for fiber breakage, but not the friction and lubrication force 
explained in Section 3.4.2.2. Although these forces can be approximated, 
their influence on fiber breakage is unclear [68,70]. The fiber breakage is 
based on buckling, and therefore on pressure. The contact force may be set in 
relation to the direction of the eigenvectors and the contact points may be 
assumed equidistant along the fiber, but still it is unclear if the forces favor or 
hinder buckling and breakage. Additionally, the Jeffery orbit is only based on 
hydrodynamic forces and this aspect would not be valid any longer. Future 
studies have to focus on the role of contact forces in fiber breakage to enable 
a meaningful use in such simulations. Furthermore, the information may be 
used for a dynamic description of the breakage coefficient Cpr, which is 
assumed to be constant to this point of time. 


The main difference of the novel approach presented within this work and the 
state-of-the-art, is the calculation of the forces and the derived buckling 
criterium. In Phelps et al. [70] the forces are calculated by the approach of 
Dinh and Armstrong [72] (see Eq. (2.37)). In this thesis, the forces are 
determined by considering hydrodynamic drag and lift force, as described in 
Section 3.4.1. These forces act on reference fibers, represented by the eigen- 
vectors of the orientation tensor. 


Within [70], -2]2;,]a:q; =1 is assumed, so the buckling criterium is 
independent of orientation and determined by matrix viscosity, shear rate, 
fiber geometry and a fitting parameter. Furthermore, the assumption 
-2|D;;]a:q; =1 leads to an overestimation of BP, since it is 
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-2|D;,]a:q; < 1 for almost all possible orientations (cf. Eq. (2.38) and 
Figure 2.12). Additionally, this approach leads to a model where all fibers 
within one cell offer potential to break. The novel approach distinguishes 
between three reference fibers in every cell, where the Jeffery orbit and 
forces are calculated independent for every reference fiber. The final break- 
ing probability is weighted with the corresponding eigenvalues. Consequent- 
ly, this novel approach offers the following two advantages. One is, that the 
amount of fibers, able to break is determined more accurate. The other one is, 
that the forces are determined without material specific fitting parameters, so 
the breakage modeling needs one empirical parameter less compared to the 
state-of-the-art. 


3.5 Modeling of warpage and residual stresses 


3.5.1 Schematic procedure 


The simulation of warpage and residual stresses, needs information of the 
mold filling simulation and material state. Therefore, pre-processing algo- 
rithms, providing and transforming the relevant data are necessary. The 
general procedure of pre-processing is schematically shown in Figure 3.6. 


The results of the mold filling simulation are mapped, including the fiber 
orientation tensor, for orientation averaging as well as the temperature and 
curing distribution, for the simulation itself. The input of material data, 
including information of the matrix in rubbery and glassy state, as well as 
mechanical and geometrical information of the fibers and the fiber volume 
content are used for the homogenization. The results of the homogenization 
are the effective mechanical and thermal properties as well as the coefficient 
vectors of thermal expansion and chemical shrinkage 9!" and 9". The exact 
methods and models, which are used for the pre-processing and the simula- 
tion, are mentioned and explained in the following sections and a summary is 
given in Table 3.2 at the end of Section 3.5.2.3. 
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Input: process simulation Input: Material data 
(A;j,T,c) (fibers, matrix rubbery & glassy) 
Homogenization 


(Ejro A1", cp, 8", Of") 
Mapped Data 


Transverse isotropic 
material model 
Orientation averaging 
Orthotropic material model 
Structural analysis 


Figure 3.6: Schematic illustration of data transformation as pre-processing of structural analysis 


Mapping 


3.5.2 Material model 


3.5.2.1 Mechanical model 


To perform an adequate simulation of residual stresses and warpage, in a 
fiber reinforced part, the material model of the solid-mechanical structural 
simulation should be anisotropic. The material model for the fiber material is 
elastic, isotropic and independent of temperature, since glass fibers are 
simulated. For modeling the mechanical behavior of the matrix, the CHILE 
model, introduced in Section 2.3.1.2 is used. The anisotropy arises due to the 
homogenization of fibers and matrix to a homogeneous FRP. Such an ap- 
proach was already successfully applied by Bernath [14] in case of continu- 
ous fiber reinforced epoxy parts, with a transversely isotropic behavior. 
Within this work, the model of Bernath is extended to represent orthotropic 
material behavior, being more meaningful for discontinuous reinforced 
injection molded parts, which are not fully aligned and hence not transversely 
isotropic. 


The mechanical behavior of the matrix is given by Eq. (2.49). The elastic 
modulus in the glassy state is assumed to be constant. In the rubbery state, the 
elastic modulus is given by 
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c2 — c2\ 8/3 
E-=E (=) ; (3.83) 
cr cr,max 1 = c 


with Eecr,max representing the elastic modulus in the fully cured rubbery state 
(c = 1) and c, being the point of gelation. This power law based approach is 
presented by Adolf and Martin [87] and also used by Bernath [14]. 


Besides the mechanical load, mainly based on holding pressure, the part 
undergoes chemical and thermal change, evoking thermal and chemical 
shrinkage. Hence, the complete change of strain is given by 


Asi = AeMe> + Aeth + Neff, (3.84) 


with change of thermal strain dein being defined by 

Ac® = HAT (3.85) 
and change of chemical strain Leit being defined by 

Acs = gf" ac, (3.86) 
where Aef? = Ae} = 0 for i + j. 


For the homogenization of the mechanical properties the method presented 
by Tandon and Weng [96] is applied to create a transversely isotropic materi- 
al, which is orientation averaged afterwards with the method presented by 
Advani and Tucker [51]. This procedure is similar to Hohberg [99]. Descrip- 
tion of oth and on are given in Section 3.5.2.2 for thermal expansion and 
3.5.2.3 for chemical shrinkage. 


3.5.2.2 Thermal model 


The thermal properties, which need to be considered for the warpage analysis 
are the specific heat capacity, thermal conductivity, coefficient of thermal 
expansion and glass transition temperature. The latter depends on the matrix 
material and, since thermoset materials are regarded, on the degree of cure. 
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The cross-linking due to curing complicates the sliding of the polymer 
molecules and hence T, rises with the degree of cure. A modeling approach 
for the relationship between T, and the degree of cure is given by DiBene- 
detto [125], by 


Ts 0 — Ts o JKToC 
Ty = Teo q (Tas — Tao)erse (3.87) 
" 1- (1 — Krg)c 
with Tyo and T,. being Ty for c = 0 and c = 1 and krg as material specific 
dimensionless modeling parameter. This approach is well established and 
applied in several studies [14,126-128]. 


The effective thermal conductivity is determined according to the work of 
Clayton [32]. Contrary to the mold filling phase, the holding pressure and 
curing phase take up a longer time period, with a significantly higher amount 
of heat flow. Therefore, it is reasonable to model the thermal conductivity as 


anisotropic vector AN. 


Within a transversely isotropic FRP, the conductivity is only distinguished in 
fiber direction Ai = Al and perpendicular A = at = a". In fiber direction 


the property is simply given by volume averaging so 
AY = Prims + (1 - Pimm (3.88) 


with Ah f and Ahm representing the thermal conductivity of fibers and 


matrix, which are assumed to be isotropic. 


For the perpendicular direction, Clayton [32] extends the general form for 
dilute dispersion given by 


ua Ang—A 
2 th,M = ( th,f th,M ) (3.89) 


AD + Wim Acne + HAtn,m 


and the corresponding differentiation 
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(u + DAnm Ath,f — AthM 
t! dath >; t t PA 


ET Ne eg (3.90) 
(at + A) Ane + Mtn 
with the dimensionless shape factor u. It is 
aay Ans \ do 
2 -= ( thf = 42 ) f (3.91) 
1+W% Arne + HAm/1- Df 


assuming AN = Ano for small additions of &, and 1/(1 — &;) as correction 
factor for this assumption. Further the integration from Ab = Ahm at Oe = 0 
to AM = Ans at Op = 1 leads to 


t 
th IFu 
Tems (a) (ae) (3.92) 


Ane + Anu) \ Ae 


The shape factor is given by u = 1 for cylindrical inclusions, such as fibers 
[32], so Eq. (3.92) leads to 


1 
th th \2 
ae (2) üs (+ _ 1)( A2 | (3.93) 
Anm Am Atı,M Ath,M 


Finally, the thermal conductivity perpendicular to the fiber direction is given 
by 


th _ 
Ay = 


A 4 ei 
EM |, aa oo ( thf _ 1) 42i 
4 Ath M 


(3.94) 


- (1-0) (7 - 1) 


Anm 
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3.5 Modeling of warpage and residual stresses 


The last property needed for warpage analysis is the thermal expansion 
coefficient. Here, the energy based approach of Schapery [97] is used. Based 
on the energy potential, bounds are defined for a linear elastic surface traction 
problem with space-wise constant temperature change. Within [97] it is 
shown, that the upper bound shows the minimal error for volume averaging, 
so 


DE Oye + (1 — Op) EMP 
DEp+(1-P)Ey 


on = (3.95) 


where Ey is the elastic modulus of the matrix. The same procedure for the 
perpendicular direction results in 


Os = OS = (14 v)O One ++ vM)Pmdınm (3.96) 

u] ; 
with vç and vy being the Poisson ratio of fibers and matrix. Similar to the 
mechanical properties, the thermal conductivity and expansion coefficient are 
only determined for a transversely isotropic material. Therefore, the proper- 
ties are further orientation averaged to create an orthotropic material with the 
orientation averaging scheme presented by Advani and Tucker [51] and 
explained in Section 2.3.2.2. 


3.5.2.3 Model for chemical shrinkage 


The homogenization of the chemical shrinkage is not implemented within 
homogenization the algorithm. Therefore, it is determined by simple volume 
averaging as 


od = G0 +1 ou, (3.97) 


with ou and gre being the coefficients of chemical shrinkage for pure 
fibers and pure matrix. Since there is no chemical shrinkage in the fibers 
9"! = 0), Eq. (3.97) simplifies to 
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9% = (1-6) ee. (3.98) 

The chemical shrinkage is assumed to be isotropic, so 
9 = (el Bl lol. (3.99) 
Table 3.2 gives a summary and overview of the homogenized mechanical and 


thermal material properties with corresponding homogenization method and 
whether it is orientation averaged or not. 


Table 3.2: Summary of homogenized material properties with corresponding homogenization 


methods 
Property Homogenization method on averag- 
ESk Tandon and Weng [96] v 
ieee Tandon and Weng [96] Vv 
a Clayton [32] Vv 
Cp Volume averaging x 
om Schapery [97] Vv 
geh Volume averaging x 
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4 Verification and validation 


4.1 Numerical verification 


4.1.1 Phase-dependent boundary condition 


To verify the boundary condition described in Section 3.2.2 a stair-like 
model, as shown in Figure 4.1, is simulated, so different amounts of material 
reach the outlets at different points of time. At the beginning, the cavity is 
filled with air and the FRP enters the model at the inlet (blue) with a constant 
velocity of U; = (8.0:10-3 0 0) m/s. Hence the cavity should be com- 
pletely filled after 5 s, if no material leaves the model. The model is meshed 
with cubic hexahedrons, having an edge length of 1 mm. 


Figure 4.1: Stair model for verification of phase-dependent boundary condition. Inlet in blue and 
three outlets with different distances to the inlet in red. Walls in transparent 


grey 
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To keep the simulation as simple as possible, the material is modeled iso- 
tropic and Newtonian with 7 = 10 Pa's and a no slip boundary condition is 
chosen for the walls. Four different formulations of the phase-dependent 
boundary condition at the outlet are compared, as shown in Figure 4.2. 


alpha FRP 
0.00 0.20 0.40 0.60 0.80 1.00 
e l | 


Interpolated switch 
(Eq. 3.4) 


ZeroGradient 


Ci 
X] 


Switchover for a > 0 Switchover for a = 1 


Figure 4.2: Simulation results after 5 s for different formulations of the phase-dependent 
boundary condition. Cut through the x,-plane of symmetry (see Figure 4.1). 
Colors visualize the value of the VOF factor a. Switchover means change from 
zeroGradient to Uzero (Eq. (3.3)) 


The case ‘ZeroGradient’ applies Eq. (3.2) all the time and independent of 
VOF factor a, so the FRP can leave the cavity without any resistance. 
‘Switchover for æ > 0’ means a spontaneous transition from zeroGradient to 
U,ero (Eq. (3.3)) as soon as æ > 0 without any intermediate state. Similar 
‘Switchover for æ = 1’ means a spontaneous transition from zeroGradient to 
Uzero as soon as a = 1. At last the interpolated formulation, as given by Eq. 
(3.4), is regarded. Figure 4.2 shows the comparison of the simulations at the 
X2-plane of symmetry after 5 s, so the cavities should be completely filled 
with FRP. The blue areas for ‘ZeroGradient’ and ‘Switchover for a = 1’ 
indicate a relative high material loss compared to the other two formulations, 
where the cavity is almost completely filled. Since all finite volumes of the 


88 


4.1 Numerical verification 


simulation mesh have identical dimensions, the results can be quantified by 
summing up all values of æ and divide the result by the total number of finite 
volumes of the domain. The results would be one for a completely filled 
cavity. The results are shown in Table 4.1. 


Table 4.1: Normalized filled volume for different formulations of the phase-dependent boundary 
condition at the outlet 


Formulation Normalized fill state 
ZeroGradient 0.8479 
Switchover for a = 1 0.9134 
Switchover for a > 0 0.9904 
Interpolated switch Eq. (3.4) 0.9913 


As expected, ‘ZeroGradient’ generates the highest material loss. The results 
of ‘Switchover for a > 0’ and Eq. (3.4) are quite similar, although the 
interpolated boundary condition is slightly better. A reason for the similarity 
is the sharp flow front with only one or two partially filled elements along the 
flow path. A more complex material behavior resulting in a torn and chaotic 
flow front would lead to more partially filled elements at the flow front and 
hence to more remaining air in the cavity for “Switchover for a > 0’. Never- 
theless, Eq. (3.4) generates the best results and is chosen for further mold 
filling simulations. 


4.1.2 Analytical verification of isotropic flow modeling 
with anisotropic viscosity tensor 


The isotropic solution of Nij, is given by an isotropic distribution of fibers, 


represented by Ay = 1/3 and A;; = 0 fori + j (see Figure 2.11), resulting in 
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mr me me 0 0 0 
mo me O 0 0 
’ iso 
min? = er ey (4.1) 
m 0 0 
sym. nise o0 
nE? 


The values for nie, nite and 7'S° are given in Table 4.2. 


Table 4.2: Values for the entries of Wier 150 (Eq. (4.1)), calculated with Eq. (3.21) 


Tensor entry Value 

nis Anıı + 24112 + 24N23 
Il 45 

mise 211 — 1212 — 1223 
JL 45 

nis? N11 + 6N12 + 6N23 
> 15 


Here, A;;xı is determined with the linear closure approximation, representing 
the exact solution for an isotropic fiber distribution [62]. For an isotropic and 
Newtonian fluid, it would be 7,2 = 23 = 7 and n11 = 37, known as the 
Trouton ratio as relation between shear and elongation viscosity, first de- 


scribed by Trouton in 1906 [129]. Assuming this, Dae can be rewritten as 


4/3 -2/3 -2/3 0 0 
4/3 -2/3 0 0 
i 4/3 0 0 
Nj = 7 i a en 
1 


sym. 


HMHooo00o0 


being the correct fourth order viscosity tensor of an isotropic and Newtonian 


90 


4.1 Numerical verification 


fluid. Nevertheless, with respect to Eqs. (3.16)-(3.18) the only two solutions 
for 712 = N23 are Pr = 0 and ®; ~ 3.383 P „ax. Since the latter is unrealistic, 
&, = 0 is chosen to compare the anisotropic tensor to an isotropic solution, 
therefore, it is 711 = 0. Of course, this is also unrealistic, since the assump- 
tion and micro-mechanical models in Section 3.3.1.2 are not valid for &; = 
0, but it is the only way to compare the approach to one of the analytical 
approaches mentioned in Section 2.2.1.2. By assuming ® = 0 it is 742 = 
N23 = Nm and 7,, = 0. Hence the scalar effective viscosity according to Eq. 


(3.26) is given by Nertiso = 4/ 5M and the complete fourth order tensor is 


given by 
4/3 —2/3 —2/3 0 0 0 
4/3 -2/3 0 0 0 
FRP,iso _ 4 4/3 0 0 0 
Nijkı = /s NM 1 0 ol (4.3) 
sym. 1 0 
1 


Hence, if implemented correctly, the simulation creates similar results to an 
analytical approach, calculated with a viscosity being 80 % of the viscosity in 
the simulations. As reference case, a Poseuille flow as described in Section 
2.2.1.2 is chosen. The parameters are given in Table 4.3. 


Table 4.3: Parameters for comparison of the viscosity tensor to a Poiseuille flow 


Parameter Value Unit 
Pressure difference 1 MPa 
Tube length 200 mm 
Tube radius 10 mm 
Viscosity 100 and 80 Pa's 


The simulation model is built up in 3D. The boundary conditions are defined 
to fit the assumptions of a Poiseuille flow. The fiber aspect ratio is set to one 
to suppress fiber re-orientation. The mesh is created with 21 hexahedral 
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elements along the diameter and 200 along the flow path. Since the simula- 
tion method is created for injection molding simulation, it is a transient 
solver. The simulation time is set to 10 s, according to the mean volume flow, 
the tube should be 50 times completely flown through at this point of time 
and the steady state is reached. 


The results are shown in Figure 4.3, represented by the velocity profile along 
the diameter. The simulation results are slightly lower than the analytical 
ones, but fit well. The maximum velocity in the center is 1.5625 m/s for the 
analytical and 1.4991 m/s for the simulation approach, giving a maximum 
deviation of about 4 %. Regarding the complete profile the mean square error 
(MSE) is 0.0026, which represents a good agreement. One reason for the 
difference is the numerical solving and discretization in the simulation. 
Additionally, the analytical approach solves only a 1D problem, while the 
simulation is still solved in full 3D. 


Analytical 
- -v- - Simulation 


Velocity in m/s 
oO 
00 


0 5 10 15 20 


Position along the diameter in mm 


Figure 4.3: Comparison of velocity profile between analytical approach (black) and simulation 
(red) in case of a Poiseuille flow 
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The analytical solution should not be seen as correct, since it also underlies 
assumptions and simplification, which are not all identical to the ones of the 
simulation. One example is the assumption of no radial flow in the analytical 
solution, which is not suppressed in the 3D-simulation model, although the 
velocity components are very low, they are not exactly zero. In summary, the 
simulation creates meaningful results, close to the analytical approach. This 
verifies the implementation in case of pure flow modeling, independent of 
anisotropic effects. 


4.1.3 Numerical verification of modeled forces and 
interaction points? 


4.1.3.1 Creation ofindividual fibers 


The approaches for calculating fiber forces, contact points and angles pre- 
sented in Sections 3.4.1 and 3.4.2 are all based on orientations of individual 
fibers and modified to apply for a homogenized material with no information 
of individual fibers. Therefore, it is meaningful to verify these approaches by 
comparing them to results, calculated with individual fibers. The second and 
fourth order orientation tensors, needed for verification, are directly comput- 
ed with individual fibers with Eq. (2.27) and Eq. (2.28). 


For verification, 22 different orientation states (OS), consisting of 500 indi- 
vidual fibers each, are considered. The OS of a single fiber is defined by the 
two angles @ and @ as illustrated in Figure 2.10. In OS 1, the possible values 
are g E [0,7] and cos(@) € [0,1] representing quasi-isotropic orientation. 
For OS with higher numbers, the window with possible values is reduced by 
rt/20 for and 1/20 for cos(@), so the fibers become more highly orientat- 
ed towards the x,-direction. Due to this rate of change, full orientation is 
given for OS 21, represented by g = 0 and cos(@) = 0. OS 22 represents 
planar isotropic orientation, so œ € [0,7] and cos(@) = 0. Although planar 


4 Similar numerical studies and results are published in [106]. 
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isotropic orientation is not typical in RIM processes, this OS is often used for 


verifications. 


The possible values of g and @ for the different orientation states are shown 
in Figure 4.4. 
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Figure 4.4: Possible values of and @ for different orientation states illustrated by colored areas 
(cf. Figure 2.10) 


Due the definition of changing the value for cos(@) in a linear way, the 
change of 0 is non-linear. By this definition the amount of fibers in x3- 
direction decreases faster, compared to the other directions x, and x. This is 
also the case in many real fiber orientations in injection molding, due to the 
thin walled part design. Figure 4.5 shows three examples of the OS 1, 10 and 
17 with corresponding orientation tensors. The single fiber orientations are 
visualized by the red dots on a unity sphere. The different fibers are created 
by using the ‘rand()’ function of Matlab R2019, creating a random value 
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between zero and one, which is multiplied with the maximum value of p and 
cos(@) in the corresponding OS. 


Orientation state 1 
*3D-quasi-isotropic’ 


173 0 0 
2 sym. 1/3 


Orientation state 10 
‘mean orientated’ 


04 0.24 0.008 
. = 0.5 -0.003 
= sym. 0.1 


Orientation state 17 


‘highly orientated’ 
0.87 0.27 0.007 
= = | 0.12 0.02 
a sym. 0.01 


Figure 4.5: Visualization of different orientation states (red dots) on a unity sphere with corre- 
sponding orientation state number and second order orientation tensor 
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4.1.3.2 Comparison of hydrodynamic forces for single fibers and 
homogenized material 


The hydrodynamic drag and lift force are calculated with Eq. (3.33), as 
hyd_f 
) 


the relation between orientation and relative velocity is given by q;. In case of 
hyd_EV 


explained in Section 3.4.1.1. For calculations with individual fibers (F; 
a homogenized material (F; ), the eigenvectors of A;,; are used as de- 
scribed in Section 3.4.1.3. Since 500 individual fibers exist for every OS, 
pet is calculated 500 times for each OS, while EN 


times, due to the three eigenvectors. To verify the novel approach, the aver- 
age force of each OS is compared, being defined as 


is calculated three 


500 


1 
ph d_av_f ) ph d_av_f 
i fii 500 A (4.4) 


n=1 


for the individual fibers. The averaged force of the homogenized approach is 
weighted with the eigenvalues so 


3 
poe = >. N (4.5) 
n=1 


The comparison of the individual and homogenized average force is shown in 
Figure 4.6, split into the drag and lift part. Two relative velocities are regard- 
ed, the first one is chosen to be AU; = 1/Y3(1 11) m/s (Figure 4.6a) and 
the second one is AU; = (1 0 0) m/s (Figure 4.6b). The material parameters 
are dp = 0.015 mm, rẹ = 100 and ny = 20 Pa's. 


The drag force components in Figure 4.6a are all identical, since the direction 
weighting of the drag force is given by AU;, whose entries are identical. The 
lift force is zero for a quasi-isotropic orientation (OS 1), because the average 
fiber direction and AU; are the same. Of course, there are lift forces, acting on 
the single and reference fibers, but they cancel each other during the summa- 
tion for the average force. This is also the fact for planar quasi-isotropic 


orientation (OS 22), where FAY z Fifa ~ 0 and only F} -3Y exists. 
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Figure 4.6: Comparison of average drag and lift force for individual fibers (f) and homogenized 
material (based on eigenvalues, EV) with different relative velocities AU; = 


1/V3 (111) m/s (a) and AU; = (1 0 0) m/s (b) 
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In detail, the force components are not exact zero, since the fibers are created 
randomly, so the orientations are not exact (planar) quasi-isotropic. 


In Figure 4.6b it is AU, = AU; = 0, hence F£"*8" = poA — g, piaga 
decreases with higher degree of orientation since the projected area of the 
fibers in the x,-plane decreases, being relevant for this AU, and is represented 
by the change of kg, as explained in Section 3.4.1.2. The lift force is zero for 
(planar) quasi-isotropic and full alignment, since these orientations are 
symmetric to AU;. Hence pee = 0 for all orientations, since the orientation 
in x3-direction is always symmetric to AU; in this case. Furthermore, it is 
Fiiftav — 0, because the lift force acts perpendicular to AU;. 

In general, the homogenized results are in exact agreement with the individu- 
al solution as shown within [106], verifying the novel approach to determine 
the average force on an arbitrary orientation state. Nevertheless, an inference 
on single fiber forces is not possible, being a disadvantage since forces on 
individual fibers may exist, but cancel each other while calculating the 
average. Such phenomena require an approach on microscopic scale and 
modeling of individual fibers. However, the average force within one cell is 
well approximated. 


4.1.3.3 Comparison of contact points for single fibers and 
homogenized material 


Three approaches are compared to the individual fiber solution in case of 
fiber-fiber contact points. One is the interaction tensor as explained in Eqs. 
(3.48) - (3.51) and [120], representing the state of the art. The other two are 
the eigenvector based (Eq. (3.54) and Eq. (3.55)) and the polynomial fit (Eq. 
(3.56)) with the values given in Table 3.1 [106]. For all approaches the 
contact points are determined with Eq. (3.47), the difference is in the calcula- 
tion of the invariants f and g of the fiber orientation distribution. For indi- 
vidual fibers Eq. (3.44) and Eq. (3.45) are used to determine f and g. The 
creation of the individual fibers and different orientation states is described in 
Section 4.1.3.1 and the second and fourth order orientation tensors are direct- 
ly computed with the individual fibers. 
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Figure 4.7: Values of the fiber orientation invariants f (a) and g (b) for the different approaches. Comparison of eigenvector approach (yel- 
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In a first step, the results of f and g are compared for the 22 different orienta- 
tion states, shown in Figure 4.7. Due to the correction factor (see Section 
3.4.2.1) the interaction tensor and eigenvector approach fit exactly, which is 
also shown analytically within [106] for specific cases. Furthermore, both 
approaches predict a too low f for higher OS, while the polynomial approach 
fits well besides OS 21. The maximum deviation of the eigenvector and 
interaction tensor approach is about 90 % at OS 20. Besides full orientation 
(OS 21), which is an unrealistic OS and will never be reached in the real 
process, the polynomial approach has a maximum deviation of about 8 % for 
OS 10. The MSE of the eigenvector and interactions tensor approach is 0.01 
and 0.001 for the polynomial approach, which is 10 times smaller. 


Regarding the approximation of g (Figure 4.7b) the interaction tensor is no 
longer considered, since only an approximation for f is given within [120]. 
Contrary to f the eigenvector method shows a maximum deviation for quasi- 
isotropic orientation of about 19 % and fits better for higher orientation. 
Again, the eigenvector approach fits perfectly for full alignment. The poly- 
nomial fit approach shows again good agreement for all OS and has maxi- 
mum deviation of about 2.4 % for planar isotropic orientation (OS 22). The 
MSE of the polynomial approach is about 0.000035 and therefore about 484 
times lower than the MSE of the eigenvector approach, being 0.015. In 
general, the polynomial fit shows the best results with the lowest MSEs and 
should also create the best results in case of calculating contact points. 


For comparison of contact points per fiber, according to Eq. (3.47), a theoret- 
ical material with & = 0.35 and the aspect ratios rp = 10 (Figure 4.8a) and 
rg = 100 (Figure 4.8b) is regarded [106]. As expected, the polynomial 
approach creates the best fitting results, with a maximum deviation of 29.4 % 
for rp = 10 and 292 % for rẹ = 100 in the unrealistic case of full alignment. 


The results of the eigenvector approach are slightly better than the interaction 
tensor for all OS, due to the consideration of g, creating an offset of exactly 
4dr(g + 1) (see Eq. (3.47)). 
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The relative importance of this offset becomes smaller for higher aspect 
ratios, since the other term to determine the contact points is linear in f and 
rg, which is the reason why 4ġ¢(g + 1) is often neglected. The maximum 
deviation of the eigenvector approach is about 29.8 % for rp = 10 and 65 % 
for rp = 100 at OS 17 and OS 19. 


The MSEs of the polynomial and eigenvector approach are 0.0917 and 1.013 
for rp = 10 and 8.83 and 82.82 for rp = 100, so the MSE of the polynomial 
approach is about 10 times smaller. In general, and as expected, the polyno- 
mial approach predicts the contact points most accurately, with respect to the 
other approaches. 


4.1.3.4 Comparison of average fiber-fiber angle and overlap area 
for single fibers and homogenized material 


In the next step, the approach to approximate the average fiber-fiber angle 
presented in [106] and Eq. (3.75) is verified. Again the 22 different orienta- 
tion states described in Section 4.1.3.1 are used. Firstly, the individual angles 
between the 500 single fibers are determined with Eq. (3.74), resulting in 
124500 angles, since one angle always includes two fibers and a fiber has no 
angles with itself (5007/2 — 500 = 124500). Afterwards, the average of 
these 124500 angles is calculated and compared to the f-based approxima- 
tion given in Eq. (3.75). The value of f is computed directly with the indi- 
vidual fibers (Eq. (3.44)) since the aim is to verify Eq. (3.75) and therefore f 
should be as exact as possible and no approximation should blur the results. 


The comparison is given in Figure 4.9, showing good agreement. Due to the 
fitting factor 4/7, there is a good agreement for (nearly) quasi-isotropic 
orientation states. The maximum deviation is about 25.8 % for OS 20 and the 
MSE is 8.45. The angle deviates about 2.25° on average. Due to the good 
agreement, the results verify Eq. (3.75) to be well suited to approximate the 
average angle within an arbitrary fiber distribution. Hence, the next step is to 
verify the modeling of the average fiber-fiber overlap area, which is based on 
this angle approximation. The comparison of the average fiber-fiber overlap 
area is given in Figure 4.10 for the aspect ratios of rp = 10 (Figure 4.10a) and 
re = 100 (Figure 4.10b). For the reference case, 124500 overlap areas, 
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according to the 124500 angles are determined with Eq. (3.73) and averaged 
for every orientation state. 
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Figure 4.9: Computed average fiber-fiber angle by averaging single angles of individual fibers 
(black, o) and approximation based on f given by Eq. (3.75) (green, V) 


Since A;j € sym., the eigenvectors are always perpendicular to each other 
and computing the individual overlap is not meaningful for the approxima- 
tions. Therefore Eq. (3.75) is used to determine the average angle and after- 
wards the average overlap with Eq. (3.73). Only the eigenvector and polyno- 
mial approach are compared, since the average angle only depends on f, the 
results of the eigenvector and interaction tensor based approaches would be 
identical. 


Both approximations predict a slightly too low overlap for both aspect ratios 
in case of lower orientation states (OS < 15). In higher orientated states, the 
eigenvector approach overshoots the results, while the polynomial approach 
fits well, except for the unrealistic state of full alignment. 
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Figure 4.10: Computed fiber-fiber overlap area for and 17 


4.1 Numerical verification 


Here, the polynomial approach shows its maximum deviation with about 
27.5 % for rp = 10 and 870 % for rp = 100, while the eigenvector approach 
fits perfectly. Again, this case will not be reached in a real process and 
therefore will not be included in the further discussion. The eigenvector 
approach has a maximum deviation of 64.5 % for rẹ = 10 (OS 17) and 162 % 
for rp = 100 (OS 20). The MSEs are 1.48-1072° for 7 = 10 and 3.36: 
1071? for rẹ = 100 in case of the polynomial approach and 6.1 102° for 
rp = 10 and 9.1-10718 for rp = 100 in case of the eigenvector approach. 
The MSE of the polynomial approach is about 4.1 times smaller for 7 = 10 
and about 27 times smaller for 77 = 100. Again, the polynomial fit approxi- 
mation creates the best results with the smallest MSE and should be used for 
further simulations, since the approximation of the lubrication force will be 
the most adequate, compared to the other mentioned approximations. 


4.1.4 Mesh study 


The interpolation of significant different material properties of FRP and air 
with the VOF in the flow front region will always create different results for 
different meshes. In general, simulations with VOF will never be independent 
of the mesh [28]. Nevertheless, the flow front shows identical tendencies, 
when reaching a corresponding mesh refinement level, such as in a mesh 
study presented in [26]. Here, the flow front of a standard isotropic Newto- 
nian and the anisotropic non-Newtonian approach are compared for different 
mesh refinements. The simulated geometry is a 3 mm thick square plate with 
100 mm edge length and a circular inlet with a 15.5 mm diameter, as shown 
in Figure 4.11. The model is meshed with hexahedral elements with 1 mm 
edge length in x,- and x-direction and 5, 10 and 15 elements in x3-direction 
(plate thickness). 


The cavity is filled with a constant volume flow of 50 cm3/s, always perpen- 
dicular to the inlet surface. Similar to [26], two different flow modeling cases 
are compared. One is isotropic and Newtonian with 7 = 100 Pa's, the other 
one is with anisotropic viscosity tensor (Eq. (3.21)) and the Castro-Macosko 
model (Eq. (2.19)) for matrix viscosity and Kamal-Malkin model (Eq. (2.11)) 
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for curing kinetics . The model parameters are given in Table 4.4 and Table 
4.5. 


Figure 4.11: CAD-model and geometry information for the mesh study. Square plate with 
circular inlet in blue and outlet in red 


The fiber orientation at the inlet is A,ı = 0.8 at the walls and A,, = 0.3 at 
the inlet center with parabolic interpolation. The further entries are A22 = 
A33 = (1-A,1)/2 and A12 = Aız = A23 = 0, if not explicitly mentioned 
differently. The fiber aspect ratio is 75. The fiber orientation is determined 
using the RSC-model (Eq. (2.30)) with x=0.03. 


Table 4.4: Model parameters for Castro-Macosko matrix viscosity model [26] 


Parameter Value Unit 
© 0.03 Pa 
NcM 0.5 - 
Com1 17 = 
er 17 - 
Bem 1.123-107 Pa's 
Tem 13750 K 

Cg 0.4 - 
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Table 4.5: Parameters for the Kamal-Malkin kinetic model [25] 


Parameter Value Unit 

R 8.3144598 J/(K-mol) 
Axmı 1.9454-10'? 1/s 

Axm2 3041.4 l/s 

Ekmı 2878805.64 J/mol 
Ekm2 38425.6452 J/mol 
MKM 1.643 - 

NKM 0.4893 - 


The simulation results of the six different cases are shown in Figure 4.12. A 
significant difference of the flow fronts between 5 and 10 elements can be 
detected in the isotropic and anisotropic case. 
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Figure 4.12: Simulated flow fronts for different simulation methods and different mesh sizes 
with increasing number of elements over plate thickness. Detail area of the 
complete geometry shown in Figure 4.11 [26] 
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Further refinement from 10 to 15 elements has almost no effect on the iso- 
tropic Newtonian results. In the anisotropic non-Newtonian case, the flow 
fronts are slightly different for 10 and 15 elements, but show the same 
tendencies, deviating from the sharp result for 5 elements. 


Both flow fronts are torn and chaotic as expected in RIM process, even 
though the mesh is structured. The torn and chaotic area is in a similar scale 
for both cases. Therefore, a mesh with 10 elements over plate thickness 
seems to reach some kind of convergence state. Although, the exact for- 
mation of the flow front will always scatter for different meshes, due to 
minimal mesh inaccuracies and numerical discretization. According to [26] a 
minimum number of 10 elements is applied between two walls in the simula- 
tions. 


4.1.5 Verification of flow and force modeling® 


In the following Sections 4.1.5.1 to 4.1.5.3 different configurations are 
simulated to highlight the benefits of anisotropic flow coupling [26] and 
show the newly gained information about the calculated forces [106]. There- 
fore, fiber aspect ratios and initial orientations are varied. The simulation 
model is identical in all cases, and similar to the model used in Section 4.1.4 
(Figure 4.11) with 10 elements in x3-direction. The material is injected with a 
constant volume flow of 50 cm?/s, perpendicular to the inlet’s surface and a 
no slip boundary condition is applied for the walls. The material has tempera- 
ture of 110 °C at the inlet and the tool temperature is 170 °C. The Castro- 
Macosko viscosity model (Eq. (2.19)) is used for the matrix viscosity and the 
Kamal-Malkin model for curing kinetics (Eq. (2.11)). The parameters are 
given in Table 4.4 and Table 4.5. The fiber aspect ratio is constant 77 = 20 or 
rg = 100, depending on the simulation. Fiber orientation model and initial 
conditions are identical to Section 4.1.4. 


5 Similar numerical studies and results published in [26] and [106]. 
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4.1.5.1 Influence of fiber length 


Figure 4.13 compares two different fiber lengths with resulting aspect ratios 
of rz = 20 and ry = 100. Due to the anisotropic viscosity tensor, depending 
on fiber length, the flow front builds up in a different way. The higher aspect 
ratio leads to higher viscosities and hence to higher viscous stress, resulting 
in a higher gradient of viscosity along the x3-direction. Consequently, the 
flow front in the simulation with longer fibers is more torn and has a larger 
area with partial wall contact. Such a flow front is a typical phenomenon of 
RIM processes as explained in Section 2.1.2.3 and would not be predicable 
within an isotropic, single-phase simulation. 


alpha FRP alpha FRP 
0.00 0.50 1.00 0.00 0.50 1.00 


t=0.2s 
t=0.3s 


rp = 20 rp=100 17 =20 rẹ = 100 


(a) (b) 


Figure 4.13: Comparison of flow fronts at the wall for two different fiber lengths: r; = 20 (left) 
and r; = 100 (right) at t = 0.2 s (a) and t = 0.3 s (b) 


4.1.5.2 Influence of fiber orientation 


Within this Section, two different initial conditions of the fiber orientation 
state are compared. The first one is the one explained in Section 4.1.4. The 
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second one is the reverse case, so it is A,, = 0.3 at the walls and A,, = 0.8 
at the inlet center. Both simulations are performed with rr = 100. 


Figure 4.14 shows the results for state of filling (Figure 4.14a) and cavity 
pressure (Figure 4.14b) after 0.5 s of filling. The different orientations have 
no major influence on the flow front evolution, but there is a significant 
influence on the in-mold pressure. Due to the high degree of orientation near 
the walls, the viscosity rises, leading to a higher pressure. In an isotropic 
standard simulation, the viscosity does not depend on the orientation state 
and therefore, no difference between the two initial states would be visible. 
This highlights a benefit of the anisotropic flow coupling, with focus on 
better pressure prediction especially near gates and other regions with signifi- 
cant change of fiber orientation. 


alpha FRP Pressure in MPa 
0.00 0.50 1.00 0.00 35.00 70.00 110.00 
we (u 


A11 = 0.8 Ay = 0.8 A11 = 0.8 A11 = 0.8 
at the wall in the core at the wall in the core 


(a) (b) 


Figure 4.14: Comparison of state of filling (a) and corresponding cavity pressure (b) at the wall 
for different initial fiber orientations. 80% fibers aligned at the wall (left) com- 
pared to 80% fibers aligned in the symmetry plane (right). The aspect ratio is 
r = 100 in both cases 
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4.1.5.3 Force distributions 


Figure 4.15 shows the A,, component results of the injection molding simu- 
lation with rẹ = 100. The resulting average hydrodynamic force (Eq. (3.33)), 
friction force (Eq. (3.66)) and lubrication force (Eq. (3.69)) are shown in 
Figure 4.16, Figure 4.17 and Figure 4.18. The average hydrodynamic force 
(Figure 4.16) is quite independent of the fiber orientation. Due to the high 
relative velocity, the force is higher near the inlet, compared to the remainder 
part. Regarding a path along the x3-direction, the force is low at the wall 


elements and in the core region, which is also traceable to the low relative 
velocity in these regions. 


Figure 4.15: Injection molding simulation result of fiber orientation tensor component A,,. Cut 
through the x,-plane of symmetry (cf. Figure 4.11) 


The comparably high forces near the wall, but not directly at the wall corre- 
spond to the regions where fiber breakage is observed in the real process. Of 
course, fiber damaging is also observed directly at the walls, but other phe- 
nomena such as fiber-wall interactions, which are not considered within this 
simulation are more important here. 
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F_hyd_avinN N 
0.00 0.02 0.4 0.06 
er 


Figure 4.16: Injection molding simulation result of average hydrodynamic force ||F,""~“"|| on a 
single fiber according to fiber orientation shown in Figure 4.15. Cut through 
the x,-plane of symmetry (cf. Figure 4.11) 


F_fr_av in N 
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Figure 4.17: Injection molding simulation result of average friction force on contact point ||” || 
according to fiber orientation shown in Figure 4.15. Cut through the x,-plane 
of symmetry (cf. Figure 4.11) 
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Contrary to hydrodynamic and lubrication forces, the friction force (Figure 
4.17) is not depending on an absolute value of any relative velocity. Hence it 
is not higher near the inlet. The friction force is depending on the specific 
number of contacts N;.g, which is higher in less orientated regions, as shown 
in Section 4.1.3.3. Therefore, the friction is higher in the core region, where 
the fibers are less orientated and lower at the walls, especially in regions with 
high degree of orientation (cf. Figure 4.15). 


The lubrication force distribution, shown in Figure 4.18, is mainly depending 
on fiber-fiber overlap and velocity gradient. Hence, the high regions are 
regions with high degree of orientation (cf. Figure 4.15), due to the larger 
overlap area. 


F_lu_avinN 
0.00e+00 5.00e-6 1.00e-05 
| OOOO A 


Figure 4.18: Injection molding simulation result of average lubrication force on contact point 
IF | according to fiber orientation shown in Figure 4.15. Cut through the x,- 
plane of symmetry (cf. Figure 4.11) 


Most high values are reached by a combination of high degree of orientation 
and high velocity gradients as detectable at walls near the inlet. Another 
region with slightly higher lubrication forces is near the walls and near the 
flow front, where the fiber orientation is almost isotropic, but the velocity 
gradient assumes high values. 
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In general, all three computed forces show meaningful results along the flow 
path with respect to fiber orientation and flow field. The friction and lubrica- 
tion force show both a significant dependency on orientation. Therefore, an 
orientation dependent modeling of these forces is essential when using them 
in further modeling approaches. Knowledge about the mentioned forces may 
help to improve state of the art models of fiber orientation, breakage, fiber- 
matrix separation or viscosity and flow modeling itself. The hydrodynamic 
forces, used to calculate fiber breakage, are experimentally validated in 
Section 4.2.3. 


4.2 Experimental validation of filling 
simulation 


4.2.1 Anisotropic flow modeling® 


The novel approach for anisotropic flow modeling should be experimentally 
validated. The validation should include fiber materials with fiber dependent 
quantifiable flow effects. Therefore the experimental work of Chiba and 
Nakamura [130] is referred. Here, the steady state flow of a Newtonian fiber 
suspension in a 1:4 backward-facing step channel is determined for different 
Reynolds numbers (Re). Caused by the backward step, eddies built up for 
higher Re, influencing the flow field and fiber orientation. The suspension is 
a Newtonian water-syrup-solvent with ny = 0.0463 Pa's, filled with ®, = 
2-10”? vinylon fibers, having an aspect ratio of rg = 283. Due to the low 
fiber volume content, the Folgar-Tucker model (Eq. (2.29)) is used for fiber 
orientation modeling. The simulation model with geometrical information is 
shown in Figure 4.19. Due to the large model, the mesh is created with cubic 
elements having an edge length of 2 mm. Hence there are areas with less than 
10 elements between two walls, which is acceptable in this case, since the 
complete cavity is filled with suspension and no flow front is regarded. 


6 Results also published in [26]. 
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Figure 4.20 represents the simulation results in case of streaming lines and 
fiber orientation ellipsoids, calculated with the eigenvalues and eigenvectors 
of the second order orientation tensor. 


Figure 4.19: Simulation model for the experimental work of [130], 1:4 backward-facing step 
channel for creating a recirculating flow. Inlet in blue, outlet in red and green 
frame, representing the area of Figure 4.20 


The results are obtained within the green frame in Figure 4.19 +2 mm in x3- 
direction, so some streamlines disappear, since they leave this control vol- 
ume. The streamlines visualize no eddy for Re=1.9, but for Re=11.1 and 
Re=14.9. The vortex increases with Re. This effect and the dimensions of the 
eddies correspond to the numerical and experimental results of Chiba and 
Nakamura. Also similar to [130], the fibers align along the streaming lines 
with higher orientation near the walls and more randomly orientation in the 
vortexes. The good agreement validates the simulation approach in case of 
flow modeling and fiber orientation for low fiber volume content. [26] 


115 


4 Verification and validation 


[97] (617 aamdıq Jo) Anourwäs zo auejd-Ex ayy ur ((9) EFT = A N TTI = 9Y *(®) 
6'1 = ey) Sloquinu spjoukoy JusısypIp 107 syInsoy '[0E1] ur Suswmedx> fesuioumu 0} s[qeredtuOd are s}[NSoyY Iosua) 
Kyısoasıa DIGO.HOSTUL YIM panus (p21) SIOSUI} UONEJUSLIO 19q17 ay) JO spIosdijjo pue (yortq) sourpueans Jo uosuedwog :Q7'p 9M3 


40 mm 


50 mm 


a) 
116 


= 
oO 


4.2 Experimental validation of filling simulation 


4.2.2 Fiber orientation’ 


4.2.2.1 Material, models and process conditions 


For validation of fiber orientation, the work of Englich [9] is considered, 
containing experimental data of a 35 %-vol short glass fiber filled phenolic 
compound. The fiber orientation is determined with microscopic images, 
presented in [9]. The images are created of RIM experiments with a square 
plate, having an edge length of 150 mm, a thickness of 2 mm and a 135 mm 
long and 1.5 mm high line inlet along one edge. The complete cavity with 
sprue is given in Figure 4.21. The sprue is meshed with 27 hexahedral ele- 
ments along the diameter and 115 in x3-direction. The plate is meshed with 
300 hexahedral elements in x,- and x-direction and 10 in x3-direction. The 
simulation is performed anisotropic, non-Newtonian and non-isothermal. The 
FRP is injected with 16 cm/s and 110°C into the mold with a constant 
surface temperature of 175 °C. Similar to [26] and Section 4.1.4, the Castro- 
Macosko viscosity model (Eq. (2.19)) is used for the matrix viscosity and the 
Kamal-Malkin model for curing kinetics (Eq. (2.11)), the parameters are 
given in Table 4.4 and Table 4.5. 


The fiber orientation is determined using the RSC-model (Eq. (2.30)) with 
strain reduction factor k=0.03 and x=0.08 for a reference case with isotropic 
viscosity modeling. Both strain reduction factors are optimized to fit the 
orientation distribution [26]. The fiber orientation at the inlet is chosen to be 
A33 = 0.8 at the walls and A33 = 0.3 at the inlet center with parabolic 
interpolation. The further entries are defined as 411 = A22 = (1 — A33)/ 
2and Aı2 = Aız = A23 = 0, being also similar to Section 4.1.4, but the 
fibers are orientated dominantly in x-direction, not x4, due to the direction 
of the sprue. The fiber aspect ratio is constant rẹ = 15 [26]. The isotropic 
reference case is also simulated with the Castro-Macosko viscosity model, 
the model parameters are given in [25]. The curing modeling of the isotropic 
case is identical to the anisotropic one, so the parameters are given in Table 
4.5. 


7 Results also published in [26]. 
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Orientation 
measurement area 


Figure 4.21: Square plate mold model for fiber orientation simulations, derived from [9]. 
Circular inlet in blue and outlet in red [26] 


4.2.2.2 Results for fiber orientation distribution 


Figure 4.22 shows the results, represented by distribution of A,, over plate 
thickness in the plate’s center. 
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Figure 4.22: Second-order orientation tensor component A,, over plate thickness for RSC-model 
with isotropic (red) and anisotropic flow modeling (green). Experimental result 
(black) derived from microscopic images in [9] [26] 
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The experimental results are given by the numerical mean with correspond- 
ing standard deviation. Both simulation approaches create good results 
validating the anisotropic approach in case of fiber orientation modeling. The 
anisotropic results fit better near walls, while the isotropic result is better in 
the core region. Near the walls (first two and last two points) the MSE is 
0.0042 for the isotropic and 0.0016 for the anisotropic simulation. In the core 
region (three middle points) the mean square error is 0.0069 for the isotropic 
and 0.019 for the anisotropic simulation. 


Regarding the complete region the mean square error is 0.0054 for the 
isotropic and 0.0091 for the anisotropic simulation. Although the results are 
quite similar, the strain reduction factor is more than 2.5 times higher in the 
isotropic simulation, highlighting the influence of fiber-flow coupling for 
fiber re-orientation. 


4.2.3 Fiber length 


4.2.3.1 Material, models and process conditions 


For the experimental investigations focusing on fiber length and cavity 
pressure within this study, an experimental long fiber phenolic molding 
compound of the novolac type is considered. The compound Porophen 
GF9201L12a by Sumitomo Bakelite Co., Ltd. was supplied as unidirectional- 
ly glass fiber reinforced flakes with an initial fiber length of Lfọ = 12 mm 
[131]. According to the measurements, the fibers are assumed with a constant 
diameter of ds = 0.0017 mm. This compound is typically used for compres- 
sion molding applications, not for injection molding, due to the processing 
difficulties. However, it is used for injection molding in a study to develop a 
long fiber RIM process in a DFG-project running since 2018, with project 
number 400343062. These experiments have been performed at the Fraunho- 
fer ICT in 76327 Pfinztal, Germany [132]. 


Processing the experimental long fiber phenolic compound has proven to be 
difficult in the injection molding process, due to the strong adhesion of the 
resin to the screw and the barrel of the injection molding machine, resulting 
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in step-like and abrupt screw movement during plasticizing. Previous rese- 
arch studies by other authors reported similar problems [133]. A higher back 
pressure and higher barrel temperatures had to be used to counteract the 
adhesion, compared to conventional short fiber molding compounds. The 
higher temperatures reduce the adhesion to the barrel, but of course also 
support the curing of the resin. In order to avoid premature curing, plasticiz- 
ing was only performed with the plasticizing unit retracted from the hot mold 
[132]. The process parameters found in this way are given in Table 4.6. A 
stable and repeatable process was possible. 


The injection molding experiments were performed on a KraussMaffei 
550/2000 GX injection molding machine equipped with a standard 60 mm 
thermoset screw and no non-return valve. The temperature control of the 
plasticizing unit is realized by four individually controlled, oil-tempered 
zones. The clamping unit has a maximum clamping force of 5500 KN. For in- 
mold pressure measurement, two capacitive sensors of the type 6163 manu- 
factured by Kistler Instrumente GmbH (Sindelfingen, Germany) are placed 
within the mold. Plates with a constant wall thickness of four millimeters and 
a size of 480 mm x 190 mm are manufactured. The mold has a central sprue 
with a cone shape and 185 mm height, further defined by a start diameter of 9 
mm, an end diameter of 15.5 mm. The simulation model of the plate with 
sprue and part of the screw chamber is shown in Figure 4.23. 


In the simulation, the front part of the plasticizing unit barrel and the nozzle 
are considered (Figure 4.23b). By doing so, a more realistic modeling of fiber 
orientation and length distribution at the beginning of the sprue is possible, 
since the initial fiber length distribution is determined in the frontal part of 
the plasticizing unit barrel (see Section 4.2.3.2). The model is meshed with 
hexahedral cells. The cells in the plate have a length of 2.3 mm in x;-, 
1.75 mm in x,- and 0.4 mm in x3-direction. The sprue and nozzle are built up 
with 24 cells along the diameter, 32 along the circumference and 1.5 mm 
edge length in x3-direction. 
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Table 4.6: Process parameters for injection molding experiments of the long fiber phenolic 
molding compound [132] 


Parameter Value Unit 
Mold temperature 185 ae 
Barrel temperature profile 50-70-90-90 °C 
Plasticizing volume 650 cm? 
Screw speed plasticizing stage 1 20 l/min 
Back pressure plasticizing stage 1 120 bar 
Plasticizing stage 1 volume 100 cm? 
Screw speed plasticizing stage 2 65 l/min 
Back pressure plasticizing stage 2 100 bar 
Plasticizing stage 2 volume 650 cm? 
Injection speed 150 cm/s 
Cavity pressure for switching to holding 100 bar 
pressure 
Holding pressure stage 1 600 bar 
Holding pressure stage 1 duration 15 s 
Holding pressure stage 2 600 — 10 (linear bar 
ramp) 
Holding pressure stage 2 duration I s 
Cure time 180 s 


The screw chamber is meshed with 72 elements along the diameter, 64 along 
the circumference and 1.57 mm edge length in x-direction. 


At the beginning of the simulation, the screw chamber is filled with FRP, 
while the rest of the cavity contains air. In the following, FRP enters the 
model at the top surface of the screw chamber, which is highlighted in blue in 
Figure 4.23. Both short edges of the plate are defined as outlet, but only one 
is visible in Figure 4.23a. The matrix viscosity, curing and fiber orientation 
model are identical to Section 4.2.2.1. Since the resin system is quite similar 
to the one considered in Section 4.2.2.1 the same parameters are used for the 
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viscosity model, but the viscosity is lowered to 50 %, since the resin system 
of the long fiber material shows a lower viscosity with same tendencies on 
temperature and shearing. 


Screw Chamber 


(b) 


Figure 4.23: Cavity for validation (rectangle plate with central sprue) with positions of pressure 
sensors p, and po, area for fiber length evaluation Ly and inlet area in blue and 
outlet in red. Complete cavity (a) and cut through x,-plane of symmetry for de- 
tailed visualization of screw chamber and nozzle (b) 
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The thermal boundary conditions are of Dirichlet type and fit to the process 
boundary conditions (see Table 4.6). For the velocity, a no-slip boundary 
condition at the walls is assumed, which is a simplification in RIM simula- 
tion, but creates good results as mentioned in Section 2.2.1.5. At the inlet the 
velocity is uMlet=(0 0 -0.425) m/s. For the velocity at the outlet, the 
phase-dependent boundary condition given by Eq. (3.4) is applied. The 
material has an initial cure of 0.01 % and the fiber orientation distribution at 
the inlet and in the initially stored material is quasi-isotropic. The fiber length 
distribution at the inlet and in the initially stored material is identical to the 
measurements in the screw chamber (see Section 4.2.3.2). All other boundary 
fields (fiber orientation at the walls, curing rate, etc.) are set to zeroGradient 
at the walls. Only the mold filling with a time of 3.3 s is simulated, so the 
mold is filled. Since no relative velocity appears in the completely filled 
mold, no fiber breakage will occur in the simulation. Therefore, ongoing 
process steps like holding and curing have no effect on the fiber length 
distribution within the simulations. 


4.2.3.2 Results for fiber length distribution 


For the fiber length analysis, specimens are extracted from the frontal part of 
screw chamber and the molded part (position L; Figure 4.23). The measure- 
ment is carried out by and according to the method described and validated 
by Maertens et al. [132]. The fiber length measurement principle is destruc- 
tive and based on the commercially available FASEP system (IDM Systems, 
Darmstadt, Germany). The measurement process works in four steps accord- 
ing to the work of Goris et al. [134]. For the first step (matrix removal), a 
sample of approximately 2.5 g is extracted from the plate and the matrix is 
removed by for 18 h in pyrolysis at 650 °C and air atmosphere using a 
TGA701 device by LECO Corporation, (St. Joseph, Michigan, USA). After- 
wards, the ash residue is transferred into an aqueous solution, which is 
diluted using an apparatus consisting of a stirrer and beaker. Therefore, a 
repeatable and controlled sample-taking out of the dilution is enabled. Con- 
trary to current state-of-the-art fiber dispersion and sample-taking methods, 
the influence of the operator on the measurement results is reduced [132]. 
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The third and fourth step are the image acquisition and the image processing 
(fiber detection) by using the algorithms provided by the FASEP software. 


For the simulation results, the fiber breakage is simulated with the method 
described in Section 3.4.3. For the child generation rate and the length 
evolution, the state-of-the-art approaches given by Eq. (2.40) and Eq. (2.44) 
are used. The force is determined with eigenvectors and hydrodynamic forces 
as described in Section 3.4.1 and based on these forces, the breakage proba- 
bility is determined with Eq. (3.79). 24 different fiber lengths between 0.25 
mm and 11.75 mm (in 0.5 mm steps) are possible, corresponding to the 
clustering of the fiber length measurement, being also 0.5 mm. The breakage 
coefficient is set to be Cpr = 0.025 and the standard deviation for breaking 
point is Spr = 1, being identical to the values used by Phelps et al. [70]. The 
difference is the force modeling, which is based on hydrodynamic forces and 
set in relation to the eigenvectors of the orientation tensor in this work (Eq. 
(3.33)). 


The experimental results of the screw chamber (black) and plate (green), as 
well as the simulation results in the plate (blue) are shown in Figure 4.24 for 
fiber lengths from 5.25 mm to 11.75 mm and Figure 4.25 for 0.25 mm to 4.75 
mm. The experimental results are the average of five measurements, illustrat- 
ed with corresponding standard deviation and results are weighted by fiber 
length. 


Only slightly more than 5 % of the initially 12 mm fibers remain at the end of 
plastification, meaning that most of the fiber breakage happens within the 
plastification and not during mold filling. There are some other longer fibers 
left with an amount of 2 % around 6 mm and 8 mm length, but most of the 
fibers (about 70 %) are already shorter than 1 mm before entering the mold. 
Regarding the plate measurement, nearly no fibers longer than 3.25 mm are 
left, since the amount is about 0.6 % in sum. Therefore, the amount of fibers 
shorter than 1 mm has risen to about 78.5 %. 


Regarding the complete measurement of screw chamber and plate, the aver- 
age fiber length in the experiments decreases by about 14 % from 0.434 mm 
to 0.38 mm. The simulation corresponds to the measurements, predicting 
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about 79.6 % to be shorter than 1 mm, being a deviation of only 1.1 %. As 
shown in Figure 4.24 and Figure 4.25 the simulated fiber length distribution 
is in good agreement with the experimental data, showing a slight rise of 
proportions from 11.25 mm to 1.75 mm and a following faster rise of propor- 
tions up to 0.25 mm. 


0.08 


Exp. screw chamber 
Exp. plate 


0.06 L X- - Sim. plate 


0.04 


0.02 


Length weighted amount of fibers 


Fiber length in mm 


Figure 4.24: Length weighted amount of fibers with lengths of 5.25 mm to 11.75 mm. Experi- 
mental results of screw chamber (black) and in plate (green) with correspond- 
ing standard deviation. Simulation results in the plate in blue x 


One exception is 11.75 mm, where the predicted amount of fibers is about 5.6 
times higher, compared to the measurements. Nevertheless, compared to the 
initial state, nearly 87 % of the 11.75 mm fibers which break due to the 
measurement are also broken within the simulation. 


Figure 4.26 shows the simulated average fiber length for the beginning and 
end of the sprue, being critical transit areas for fiber breakage. 
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Figure 4.25: Length weighted amount of fibers with lengths of 0.25 mm to 4.75 mm. Experi- 
mental results of screw chamber (black) and in plate (green) with correspond- 
ing standard deviation. Simulation results in blue x 


By transition from the screw chamber to the sprue (Figure 4.26a), the average 
fiber length decreases already about 5 %. A state-of-the-art simulation often 
starts at the beginning of the sprue and the screw chamber and nozzle are 
neglected. As can be seen, a proportion of fibers is already broken at this 
point and the assumption of a fiber length distribution similar to the stored 
material is not valid at this point. 


Figure 4.26b further highlights the transit from sprue to mold as an important 
area of fiber breakage, where the average fiber length is reduced by 10 %, 
due to the sharp edge, creating high shearing rates and relative velocities. 
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Av. fiber length in mm 
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0.407 
0.382 
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b) 0.334 


Figure 4.26: Simulated average fiber length at the end of filling, weighted by number. Image 
detail of transition from screw chamber to sprue (a) and sprue to mold plate (b) 


The simulation results are in good agreement with the experimental data, 
validating the novel approach of simulating fiber breakage based on eigen- 
vectors and hydrodynamic forces on macroscopic scale. The fiber breakage is 
predicted at meaningful areas within the part. Of course, shortening of the 
fibers affects the anisotropic viscosity modeling, depending on the fiber 
aspect ratio. Hence, the simulated cavity pressure is also influenced by the 
fiber breakage, which will be discussed in the next Section. 


4.2.4 Cavity pressure 


Figure 4.27 shows the experimental and simulation results of cavity pressure 
for two positions within the mold, as illustrated in Figure 4.23. The experi- 
mental pressure data is again the average of five measurements with corre- 
sponding standard deviation. The data fits to the fiber length measurements 
presented in Section 4.2.3.2 so the process parameters are given by Table 4.6. 
Two simulations are compared. Both are based on the anisotropic viscosity 
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tensor presented in Section 3.3. One simulation is calculated with fiber 
breakage, where the fiber length distribution is given in Section 4.2.3.2. 
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1- - - -Sim. break p2 


Pressure in MPa 
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Figure 4.27: Cavity pressure during mold filling at sensor p, (solid lines) and p, (dotted lines) as 
shown in Figure 4.23. Comparison of experiments (black) and anisotropic 
modeling without fiber break (red) and with fiber break (green). Experimental 
lines are the average of five measurements with corresponding standard devia- 
tion 


The other simulation assumes a constant fiber length of 0.434 mm, based on 
the measurements in the screw chamber. Since the simulation model is 
identical to the fiber length simulation, all parameters and boundary condi- 
tions are given in Section 4.2.3.1. Both simulations predict the cavity pres- 
sure well at position pı and slightly too high at position p, being still within 
the standard deviation of the experiments. The results match the previous 
work [26], where the anisotropic flow coupling with constant fiber length is 
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validated to be suitable for in-mold pressure simulation during mold filling. 
The focus point of this comparison is the deviation at p, of the two simula- 
tions after 4.5 s, where the simulation without fiber breakage predicts a 
higher pressure. The maximum deviation is about 0.43 MPa at 5.45 s, being 
about 10 %. The deviation is within an area, where simulated pressure is too 
high, compared to the experimental results, so in summary, the results with 
fiber breakage are slightly better than without. 


As expected, the predicted pressure, or pressure drop between pı and p, is 
higher without fiber breakage, due to the longer fibers, causing a higher 
viscosity. Although the difference is quite marginal with about 10 %, it 
corresponds to only a slight change of fiber length with about 14 % (see 
Section 4.2.3.2) and highlights the benefit of a fiber length dependent viscosi- 
ty modeling with parallel fiber breakage calculation. However, the predicted 
pressure at p, is too high for both simulations. One reason may be the no-slip 
boundary condition, which should be further investigated. Sliding along the 
wall, would lower the pressure by lowering the kinetic energy, needed for 
material transport. Of course, this would affect both regarded positions pı 
and p,, but pz will be affected stronger, due to the longer flow path. 


4.3 Outlook on prediction of warpage 
simulation 


This Section will give a short outlook on the potential of predicting warpage 
and residual stresses, presented in Section 3.5. Due to a lack of information 
about material data and experimental data, the approach cannot be validated 
within this work. But at least the general feasibility of such a simulation will 
be shown. 
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4.3.1 Modeland procedure 


4.3.1.1 Initial state of warpage analysis 


The warpage analysis begins after mold filling, including the process steps 
holding, curing, ejection and out-of-mold cooling. Four different cases are 
considered. On the one side, the time period of the curing step is performed 
for 10 s and 60 s, on the other side both configurations are simulated with an 
ejection step and without. The warpage analysis is performed with the FEM 
software Abaqus 2018 (Dassault Systémes, Johnston, RI, USA). 


The simulation model is a rectangular plate, similar to the one shown in 
Figure 4.23, except the thickness of the plate is only 3 mm not 4 mm so the 
model is similar to [26], where the anisotropic flow modeling without fiber 
breakage is validated. Here, fiber breakage is ignored, since the algorithm for 
homogenization and material creation can only handle constant fiber aspect 
ratios. Therefore, the short fiber material and mold used in [26] are regarded. 
As input for the warpage analysis, fiber orientation, temperature and degree 
of curing distribution at the end of the mold filling simulation are mapped, 
using MPCCI MapLib [100,101]. Therefore, all relevant parameters for 
material modeling are identical to the mold filling simulation in the initial 
state. 


4.3.1.2 Material modeling and boundary conditions 


During the holding and curing step, displacement is blocked on the whole 
part’s surface, since it is still in mold in the real process. The surface temper- 
ature is set constant to 170 °C. The holding stage applies for 40 s, the curing 
step for 10 s or 60s. Afterwards, an ejection step of 3 s is performed, fol- 
lowed by an out-of-mold cooling step of 6000 s. During the latter two steps, 
only the top end of the sprue is blocked for displacement to fix the model for 
numerical stability. Over the course of the ejection step, the corner nodes of 
the plate are displaced 20 mm in x3-direction during the complete ejection 
step. In the cooling step, the corner nodes have no specific displacement 
boundary condition, similar to the rest of the model. During ejection and 
cooling, a convection boundary condition for temperature is applied on the 
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whole part. The corresponding film coefficient and sink temperature are 
15 W/(m?K) and 20 °C. The effective material properties are determined with 
the schemes listed in Table 3.2 and explained in Section 2.3.2 and Section 
3.5.2. For the mechanical model, the CHILE model (Section 2.3.1.2) is used 
in an orthotropic formulation. The model parameters are given in Table 4.7 
and Table 4.8. The part is clustered in 30 different material sections, depend- 
ing on orientation state, so similar orientation states are summarized. The 
sections must not be contiguous, but have similar fiber orientations, which 
are normally nearby. According to the 30 material sections, 30 orthotropic 
materials are created, describing the material behavior also for the corre- 
sponding fiber orientation state. This procedure is performed to lower the 
numerical effort and increase numerical stability. 


Table 4.7: Material parameters of fiber and matrix for material modeling of warpage analysis 


Parameter Value Unit 

PM 1270 kg/m? 

Pr 2500 kg/m? 
Eu,g 3500 MPa 
Eur 233.33 MPa 

Er 8-104 MPa 

VM 0.35 - 

Ve 0.22 - 

DinMg 104 m/K 
Üth,M,r 3.3-10° m/K 
Inf 5.4:10° m/K 

Uch M 2:104 - 

Dene 2:104 - 

AthM 0.38 W/mK) 
Ath f 0.15 W/mK) 
Cp M 1850 J/(kg-K) 
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Cag 1850 J/(kg-K) 


Due to a lack of material data, the parameters are not correct for the rubbery 
state and glass transition of the matrix for the corresponding material, but 
orientated at the work of Bernath [14]. Nevertheless, the aim is a proof of 
concept, not generating quantifiable results. The values for chemical shrink- 
age and specific heat capacity are known for the compound and therefore 
chosen identical for matrix and fiber for simplification, to fit to the value of 
the compound. For the homogenization, a fiber volume content of 0.35 and 
an aspect ratio of 15 is applied. 


Table 4.8: Glass transition parameters of the matrix for material modeling of warpage analysis 


Parameter Value Unit 
Ta P K 
Te 3 K 
T,o 74 IC 
Tg, 224 °C 
Kg 0.268 - 


4.3.2 Results of the warpage analysis 


The results of the warpage analysis are shown in Figure 4.28 for 60 s curing 
time and Figure 4.29 for 10 s curing time. In both figures, the deformation on 
the part is scaled up five times, so the general shape of the final part is 
detectable. After 60 s of curing, the degree of cure is above 0.75 at every 
point of the part and above 0.9 within the plate. Therefore, it is T, > 190 °C 


within the plate and the behavior is linear elastic at the point of ejection. 
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a) Displacement_mag in m 
j +7.45e-03 


b) Displacement_mag in m 
i +7.45e-03 


Figure 4.28: Magnitude of displacement after 40 s holding, 60 s curing, ejection and 6000 s 
cooling. Simulation without ejection deformation (a) and with ejection 
deformation (b) 
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Figure 4.29: Magnitude of displacement after 40 s holding, 10 s curing, ejection and 6000 s 
cooling. Simulation without ejection deformation (a) and with ejection 
deformation (b) 
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Due to the linear elastic behavior for the complete part, the simulation results 
with ejection deformation (Figure 4.28b) are identical to the results without 
ejection deformation (Figure 4.28a), because the ejection deformation simply 
goes back to its initial state during cooling. Therefore, the ejection does not 
influence the warpage in this case and the deformations must have a different 
source. The displacement shown within Figure 4.28 is caused by frozen-in 
residual stresses during curing and thermal shrinkage in combination with 
anisotropic behavior. The maximum displacement is about 7.5 mm at the 
edge with maximum x} -value and a convex shell is built. 


In case of only 10 s curing time, the results with ejection deformation (Figure 
4.29b) and without (Figure 4.29a) are different, since the material does not 
behave linear elastic during ejection. When ejecting the part after only 10 s 
curing itis T, < 170 °C in most regions of the part, therefore the behavior is 
entropy elastic (rubbery state). Hence, a part of the ejection deformation is 
captured due to the change of material behavior during ejection and the final 
deformation is higher in Figure 4.29b, compared to Figure 4.29a. 


Comparison of Figure 4.28a and Figure 4.29a, shows lower deformations for 
lower curing times, resulting in a maximum displacement of only about 5.5 
mm. This is caused by less frozen-in stresses at the point of ejection, since 
the material is in rubbery state and no residual stresses can build up. The 
general distribution of deformation and final part shape is quite similar in 
both cases, since is mainly driven by thermal shrinkage and anisotropy, 
where the latter is identical in all simulations. In reality it would be unlikely, 
that a part with such small curing time, and hence degree of cure will deform 
less. Effects such as movement of the part and gravity will deform the still 
visco-elastic and not completely form-stable part, but are neglected within the 
simulation. 


In summary, the orthotropic CHILE model in combination with mapping of 
process fields is able to capture warpage of fiber reinforced injection molding 
parts, depending on process parameters, fiber orientation distribution and 
curing reaction. The displacements are in a realistic scale. Nevertheless, the 
material parameters in the area of the rubbery state and glass transition are 
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only estimated and not based on experimental data. Further investigations on 
material parameters and warpage experiments with injection molding parts 
must be performed to validate the method. 
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5.1 Conclusion 


Novel simulation approaches for fiber reinforced injection molding have been 
presented. The state-of-the-art is extended by a multiphase approach, consid- 
ering air and FRP within the simulation. A phase-dependent boundary condi- 
tion enables that the air may leave the mold, while the FRP is blocked. The 
boundary condition is formulated to minimize the material loss during simu- 
lation. 


To consider fiber length and orientation within the flow simulation, the 
viscosity is modeled non-Newtonian with a fourth order viscosity tensor, 
taking fiber orientation, length and volume fraction into account. Numerical 
studies show the influence of fiber length and orientation on flow front 
evolution and on in-mold pressure, highlighting the benefits of the tensorial 
approach. 


New approaches approximate hydrodynamic forces, fiber-fiber contact points 
and fiber-fiber forces within the homogenized material with information 
provided by the second order orientation tensor. The contact point and force 
modeling is verified with numerical experiments, including randomly gener- 
ated individual fibers. The novel approaches show good results for different 
orientations states and fiber lengths. 


The hydrodynamic forces are used to model fiber breakage with orientation 
and flow dependency during mold filling. The breakage modeling is validat- 
ed with fiber length distributions and pressure data of injection molding 
experiments. The numerical results are in good agreement with the experi- 
mental data. In combination with the viscosity tensor, the fiber breakage 
directly influences the flow modeling, which is also detectable for the simu- 
lated in-mold pressure. 
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After mold filling, relevant process data is mapped to perform a non- 
isothermal structural analysis with respect to fiber orientation, length and 
curing distribution. Therefore, the part is clustered into a defined number of 
different materials, which are all orthotropic and the values are determined by 
homogenization and orientation averaging. The simulation results are mean- 
ingful and show a prediction of warpage with respect to non-isotropic materi- 
al behavior, curing and glass transition in a realistic scale. 


In summary, the novel approaches improve the quality and level of detail for 
injection molding simulations with fiber reinforced polymers. A more de- 
tailed prediction of air traps, flow front evolution, cavity pressure and materi- 
al internal forces, depending on fiber orientation and length distribution and 
fiber volume fraction is enabled. In combination with the subsequent warp- 
age prediction, these approaches help to improve the final part’s quality. 
Furthermore, the higher degree of knowledge during the form filling supports 
the process and product design in an early stage of development, which 
decreases the loops of part design optimization and tool revision and there- 
fore has a positive effect on costs, economy and environment. 


5.2 Outlook 


Although the presented novel approaches improve the quality of an injection 
molding simulation, there is still unused potential for more improvement 
which needs to be regarded and validated. This concerns all process steps 
mentioned within this work. 


The calculated fiber-fiber interaction points and forces have no effect on the 
simulation at this point of time, which is of course not the case in the real 
process. The mentioned aspects may be used to improve fiber breakage 
calculations by more detailed force modeling or dynamic description of 
model parameters, used within the breakage model. Furthermore, they offer 
potential for a more detailed fiber orientation modeling and even flow model- 
ing itself. This applies for viscosity modeling and for the momentum balance 
equation itself, since the interaction forces represent material internal forces. 
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Nevertheless, the qualitive and quantitative influence of fiber-fiber contact is 
not well studied at this point of time and the resulting consequences on fiber 
breakage and material behavior are not completely known. Therefore, further 
investigations from experimental and numerical side are needed to improve 
the state of knowledge on these phenomena. Of course, these effects are on 
microscopic scale and cannot be separated from each other, creating a high 
experimental effort. 


Additionally, some boundary conditions can be improved, by a more detailed 
approach. This concern the wall-slip effect, which should be further regarded 
experimental and simulative. Influence of tool, temperature, pressure, vis- 
cosity and of course velocity and material on this friction contact should be 
quantified and taken into account in the simulations. Another aspect is the 
temperature boundary condition, which is often chosen homogeneous and 
constant at state of the art, or the tool with cooling/heating channels is also 
part of the simulation model, creating higher numerical effort. A more 
detailed description of heat flow as boundary condition improves the temper- 
ature modeling without significant rise of numerical effort, since no tool must 
be regarded. Here also, the heat transfer coefficient should be modeled with 
respect temperature, pressure, curing state, etc. Both extensions of the bound- 
ary conditions may improve the simulation results, especially in case of in- 
mold pressure and temperature modeling. 


Furthermore, more investigation on the early state of warpage analysis shown 
within this work should be done. In a first step the shown approaches must be 
validated (or modified) based on experimental data, which also includes 
experiments to determine relevant material properties. Furthermore, the 
method can be improved by more detailed material modeling. The fiber 
aspect ratio should be included into the mapping and homogenization pro- 
cess, so the warpage analysis can consider the length distribution. 


139 


6 


[10 


= 


Literature 


Osswald TA, Beaumont JP. Injection molding handbook. 2nd ed. 
Munich: Hanser; 2008. 

Höer M. Einfluss der Material- und Verarbeitungseigenschaften von 
Phenolharzformmassen auf die Qualität spritzgegossener Bauteile. 
Chemnitz: Universitätsverlag Chennitz; 2014. 

Domininghaus H, Elsner P, Eyerer P, Hirth T. Kunststoffe. Berlin, 
Heidelberg: Springer Berlin Heidelberg; 2012. 

Eyerer P, Elsner P, Hirth T. Polymer Engineering: Technologien und 
Praxis. Berlin, Heidelberg: Springer; 2008. 

Henning F, Moeller E. Handbuch Leichtbau: Methoden, Werkstoffe, 
Fertigung. Miinchen: Hanser; 2011. 

Victrex plc. Datasheet Victrex PEEK 450CA30. [June 08, 2020]; 
Available from: https://www.victrex.com/- 
/media/downloads/datasheets/victrex_tds_450ca30.pdf?rev=33e4bc128 
a7d47e997894a0fe973cb61. 

Keyser Hd, Berg L-F, Stammler J. Composite materials ready to replace 
aluminium in internal combustion engines. JEC composites 
2016(53):32-5. 

Langheck A, Reuter S, Saburow O, Maertens R, Wittemann F, Berg LF 
et al. Evaluation of an Integral Injection Molded Housing for High 
Power Density Synchronous Machines with Concentrated Single-Tooth 
Winding. 2018 8th International Electric Drives Production Conference 
(EDPC): IEEE, Schweinfurt (Germany) 2018:1-6. 

Englich S. Strukturbildung bei der Verarbeitung von glasfasergefüllten 
Phenolformaldehydharzformmassen. Chemnitz: Universitätsverlag 
Chennitz; 2015. 

Ohta T., Yokoi H. Visual analysis of cavity filling and packing process 
in injection molding of thermoset phenolic resin by the gate- 
magnetization method. Polym. Eng. Sci. 2001;2001(41):806-19. 


141 


6 Literature 


[11] 


[12] 


[13] 


[16] 


[17] 


[18] 


142 


Tran NT, Englich S, Gehde M. Visualization of wall slip during ther- 
moset phenolic resin injection molding. Int J Adv Manuf Technol 
2018;95(9-12):4023-9. 

Michaeli W., Hunold D., Kloubert T. Formfüllung beim Spritzgießen 
von Duroplasten. Plastverarbeiter 1992;43(12):42-6. 

Kurt M, Saban Kamber O, Kaynak Y, Atakok G, Girit O. Experimental 
investigation of plastic injection molding: Assessment of the effects of 
cavity pressure and mold temperature on the quality of the final prod- 
ucts. Materials & Design 2009;30(8):32 17-24. 

Bernath A. Numerical prediction of curing and process-induced distor- 
tion of composite structures: Doctoral Thesis, Karlsruher Institut fiir 
Technologie (KIT); 2020. 

Sanchez R, Aisa J, Martinez A, Mercado D. On the relationship be- 
tween cooling setup and warpage in injection molding. Measurement 
2012;45(5):105 1-6. 

Huang M-C, Tai C-C. The effective factors in the warpage problem of 
an injection-molded part with a thin shell feature. Journal of Materials 
Processing Technology 2001;110(1):1-9. 

Choi D-S, Im Y-T. Prediction of shrinkage and warpage in considera- 
tion of residual stress in integrated simulation of injection molding. 
Composite Structures 1999;47(1-4):655-65. 

Gao Y, Wang X. Surrogate-based process optimization for reducing 
warpage in injection molding. Journal of Materials Processing Technol- 
ogy 2009;209(3):1302-9. 

Ozcelik B, Erzurumlu T. Comparison of the warpage optimization in 
the plastic injection molding using ANOVA, neural network model and 
genetic algorithm. Journal of Materials Processing Technology 

2006; 171(3):437-45. 

Batchelor GK. An introduction to fluid dynamics. Cambridge: Cam- 
bridge University Press; 2000. 

Autodesk. Autodesk Simulation Moldlfow Insight: Fundamentals 2014. 
Student Guide - Theory and Concepts Manual. Charlottesville, VA, 
USA: Ascent; 2013. 

Autodesk. Autodesk Simulation Moldlfow Insight: Advanced Flow 
2014. Student Guide - Therory and Concepts Manual. Charlottesville, 
VA, USA: Ascent; 2013. 


6 Literature 


[23] 


Tamil J, Ore SH, Gan KY, Ore SH, Gan KY, Bo YY et al. Molding 
Flow Modeling and Experimental Study on Void Control for Flip Chip 
Package Panel Molding with Molded Underfill Technology. Interna- 
tional Symposium on Microelectronics 2012 // 2011;2011(1):673-82. 
Ospald F. Numerical Simulation of Injection Molding using Open- 
FOAM. Proc. Appl. Math. Mech. 2014;14(1):673—4. 

Wittemann F, Maertens R, Bernath A, Hohberg M, Kärger L, Henning 
F. Simulation of Reinforced Reactive Injection Molding with the Finite 
Volume Method. J. Compos. Sci. 2018;2(1):5. 

Wittemann F, Maertens R, Karger L, Henning F. Injection molding 
simulation of short fiber reinforced thermosets with anisotropic and 
non-Newtonian flow behavior. Composites Part A: Applied Science and 
Manufacturing 2019;124:105476. 

Meyer N, Saburow O, Hohberg M, Hrymak AN, Henning F, Karger L. 
Parameter Identification of Fiber Orientation Models Based on Direct 
Fiber Simulation with Smoothed Particle Hydrodynamics. J. Compos. 
Sci. 2020;4(2):77. 

Deshpande SS, Anumolu L, Trujillo MF. Evaluating the performance of 
the two-phase flow solver interFoam. Comput. Sci. Disc. 2012;5(1):14— 
6. 

Kamal MR, Ryan ME. The behavior of thermosetting compounds in 
injection molding cavities. Polym. Eng. Sci. 1980;20(13):859-67. 
Kugele D. Experimentelle und numerische Untersuchung des Abkühl- 
verhaltens thermoplastischerGelegelaminate in der Prozesskette. Doc- 
toral Thesis. Karlsruhe, Germany: Karlsruher Institute of Technology 
(KIT); 2020. 

Mottram JT, Taylor R. Thermal Conductivity of Fibre-Phenolic Resin 
Composites. Part II: Numerical Evaluation. Composites Science and 
Technology 1987;29:211-32. 

Clayton WA. Constituent and composite thermal conductivities of 
phenolic-carbon and phenolic-graphite ablators. AIAA/ASME 12th 
Structures, Strucutal Dynamics and Materials Conference, Anahheim 
(California) 1971(71):380. 

Kamal MR, Sourour S. Kinetics and thermal characterization of ther- 
moset cure. Polym. Eng. Sci. 1973;13(1):59-64. 


143 


6 Literature 


[34] 


[35] 


[36] 


[37] 


[42] 


[43] 


[44] 


144 


Bernath A, Kärger L, Henning F. Accurate Cure Modeling for Isother- 
mal Processing of Fast Curing Epoxy Resins. Polymers 2016;8(11):390. 
Karkanas PI, Partridge IK. Modelling the cure of a commercial epoxy 
resin for applications in resin transfer moulding. Polymer International 
1996;41:183-91. 

Ruiz E, Waffo F, Owens J, Billotte C, Trochu F. Modeling of resing 
cure kinetics for molding cycle optimization. The 8th International 
Conference on Flow Processes in Composite Materials (FPCM8), Dou- 
ai (France) 2006:221-30. 

Meyer N, Schöttl L, Bretz L, Hrymak AN, Kärger L. Direct Bundle 
Simulation approach for the compression molding process of Sheet 
Molding Compound. Composites Part A: Applied Science and Manu- 
facturing 2020;132:105809. 

Sommer DE, Favaloro AJ, Pipes RB. Coupling anisotropic viscosity 
and fiber orientation in applications to squeeze flow. Journal of Rheo- 
logy 2018;62(3):669-79. 

Ostwald W. Ueber die Geschwindigkeitsfunktion der Viskositt dis- 
perser Systeme. I. Kolloid-Zeitschrift 1925;36(3):157-67. 

Herschel WH, Bulkley R. Konsistenzmessungen von Gummi- 
Benzollösungen. Kolloid-Zeitschrift 1926;39(4):291-300. 

Cross MM. Rheology of non-Newtonian fluids: A new flow equation 
for pseudoplastic systems. Journal of Colloid Science 1965;20(5):417— 
37. 

Castro JM, Macosko CW. Studies of mold filling and curing in the 
reaction injection molding process. American Institute of Chemical En- 
gineers Journals 1982;28(2):250-60. 

Gibson AG. Die entry flow of reinforced polymers. Composites 
1989;20(1):57-64. 

Rodgers PA. Pressure-volume-temperature relationships for polymeric 
liquids: A review of equations of state and their characteristic parame- 
ters for 56 polymers. J. Appl. Polym. Sci. 1993;48(6):1061-80. 
Padilha Jünior EJ, Soares RdP, Cardozo NSM. Analysis of equations of 
state for polymers. Polímeros 2015;25(3):277-88. 

Wang J, Hopmann C, Schmitz M, Hohlweck T, Wipperfürth J. Model- 
ing of pvT behavior of semi-crystalline polymer based on the two- 


6 Literature 


domain Tait equation of state for injection molding. Materials & Design 
2019;183(1):108149. 

Jeffery GB. The motion of ellipsoidal particles immersed in a viscous 
fluid. Proc. R. Soc. London 1922; 102:161-79. 

Advani SG. Flow and rheology in polymer composites manufacturing. 
Amsterdam, New York: Elsevier; 1994. 

Petrie JS. The rheology of fibre suspensions. Journal of Non-Newtonian 
Fluid Mechanics 1999;87:369-402. 

Folgar F., Tucker C.L. Orientation Behavior of Fibers in Concentrated 
Suspensions. Journal of Reinforced Plastics and Composites, vol. 3 
1984:98-119. 

Advani SG, Tucker CL. The Use of Tensors to Describe and Predict 
Fiber Orientation in Short Fiber Composites. Journal of Rheology 
1987;31(8):75 1-84. 

Wang J, O’Gara JF, Tucker CL. An objective model for slow orienta- 
tion kinetics in concentrated fiber suspensions: Theory and rheological 
evidence. Journal of Rheology 2008;52(5):1179-200. 

Phelps JH, Tucker CL. An anisotropic rotary diffusion model for fiber 
orientation in short- and long-fiber thermoplastics. Journal of Non- 
Newtonian Fluid Mechanics 2009;156(3):165—76. 

Skrabala O. BC. Orientierungsverhalten plättchenförmiger Zusatzstoffe 
in Kunststoffsuspensionen. Journal of Plastic Technology:157—83. 
Tseng H-C, Chang R-Y, Hsu C-H. Phenomenological improvements to 
predictive models of fiber orientation in concentrated suspensions. 
Journal of Rheology 2013;57(6):1597-631. 

Tseng H-C, Chang R-Y, Hsu C-H. The use of principal spatial tensor to 
predict anisotropic fiber orientation in concentrated fiber suspensions. 
Journal of Rheology 2018;62(1):313-20. 

Ranganathan S, Advani SG. Fiber—fiber interactions in homogeneous 
flows of nondilute suspensions. Journal of Rheology 1991;35(8):1499— 
922; 

Folgar F. Fiber orientation in concentrated suspensions: A predictive 
model. Phd. thesis. Urbana-Campaign: University of Illinois; 1983. 
Bay R. Fiber Orientation in Injection Molded Composites: A Compari- 
son of Theory and Experiment. Phd. thesis. Urbana-Campaign: Univer- 
sity of Illinois; 1991. 


145 


6 Literature 


[60] 


[61] 


146 


Phan-Thien N, Fan X-J, Tanner RI, Zheng R. Folgar—Tucker constant 
for a fibre suspension in a Newtonian fluid. Journal of Non-Newtonian 
Fluid Mechanics 2002;103(2-3):25 1-60. 

Heinen K. Mikrostrukturelle Orientierungszustände strömender Poly- 
merlösungen und Fasersupensionen. Doctoral Thesis. Dortmund, Ger- 
many: Universtität Dortmund; 2007. 

Advani SG, Tucker CL. Closure approximations for three-dimensional 
structure tensors. Journal of Rheology 1990;34(3):367—86. 

Chung, Kwon. Fiber orientation in the processing of polymer compo- 
sites. Korea-Australia Rheology Journal 2002;14(4):175-88. 

Du Chung H, Kwon TH. Invariant-based optimal fitting closure approx- 
imation for the numerical prediction of flow-induced fiber orientation. 
Journal of Rheology 2002;46(1):169—94. 

Shon K, Liu D, White JL. Experimental Studies and Modeling of 
Development of Dispersion and Fiber Damage in Continuous Com- 
pounding. International Polymer Processing 2005;20(3):322-31. 
Inceoglu F, Ville J, Ghamri N, Pradel JL, Durin A, Valette R et al. 
Correlation between processing conditions and fiber breakage during 
compounding of glass fiber-reinforced polyamide. Polym. Compos. 
2011;32(11):1842-50. 

Turkovich R von, Erwin L. Fiber fracture in reinforced thermoplastic 
processing. Polym. Eng. Sci. 1983;23(13):743-9. 

Hernandez JP, Raush. T., Rios A, Strauss S, Osswald TA. Theoretical 
analysis of fiber motion and loads during flow. Polym. Compos. 
2004;25(1). 

Durin A, Micheli P de, Ville J, Inceoglu F, Valette R, Vergnes B. A 
matricial approach of fibre breakage in twin-screw extrusion of glass 
fibres reinforced thermoplastics. Composites Part A: Applied Science 
and Manufacturing 2013;48:47-56. 

Phelps JH, Abd El-Rahman AI, Kunc V, Tucker CL. A model for fiber 
length attrition in injection-molded long-fiber composites. Composites 
Part A: Applied Science and Manufacturing 2013;51:11-21. 

Burgers JM. On the Motion of Small Particles of Elongated Form. 
Suspended in a Viscous Liquid. In: Nieuwstadt FTM, Steketee JA, edi- 
tors. Selected Papers of J. M. Burgers. Dordrecht: Springer Nether- 
lands; 1995, p. 209-280. 


6 Literature 


[72] 


[73] 


Dinh SM, Armstrong RC. A Rheological Equation of State for Semi- 
concentrated Fiber Suspensions. Journal of Rheology 1984;28(3):207- 
27. 

Titomanlio G, Jansen KMB. In-mold shrinkage and stress prediction in 
injection molding. Polym. Eng. Sci. 1996;36(15):2041-9. 

Zoetelief WF, Douven LFA, Housz AJI. Residual thermal stresses in 
injection molded products. Polym. Eng. Sci. 1996;36(14):1886—96. 
Bushko WC, Stokes VK. Solidification of thermoviscoelastic melts. 
Part I: Formulation of model problem. Polym. Eng. Sci. 

1995;35(4):35 1-64. 

Bushko WC, Stokes VK. Solidification of thermoviscoelastic melts. 
Part II: Effects of processing conditions on shrinkage and residual 
stresses. Polym. Eng. Sci. 1995;35(4):365—83. 

Kabanemi KK, Crochet MJ. Thermoviscoelastic Calculation of Residu- 
al Stresses and Residual Shapes of Injection Molded Parts. International 
Polymer Processing 1992;7(1):60—70. 

Kabanemi KK, Vaillancourt H, Wang H, Salloum G. Residual stresses, 
shrinkage, and warpage of complex injection molded products: Numer- 
ical simulation and experimental validation. Polym. Eng. Sci. 
1998;38(1):21-37. 

Bergstrom J. Mechanics of solid polymers: Theory and computational 
modeling. Amsterdam: William Andrew is an imprint of Elsevier; 2015. 
Mirabedini AS, Karrabi M, Ghasemi I. Viscoelastic behavior of 
NBR/phenolic compounds. Iran Polym J 2013;22(1):25-32. 

Prony R. Essai Expérimental et Analytique: sur les lois de la dila- 
tabilité de fluides élastique et sur celles de la force expansive de la 
vapeur de l'alkool, à différentes températures. Journal de l'École Poly- 
technique Floréal et Plairial 1795; 1(22):24-76. 

Liebig WV, Jackstadt A, Sessner V, Weidenmann KA, Karger L. 
Frequency domain modelling of transversely isotropic viscoelastic fi- 
bre-reinforced plastics. Composites Science and Technology 
2019;180(3):101-10. 

O'Brien DJ, Mather PT, White SR. Viscoelastic Properties of an Epoxy 
Resin during Cure. Journal of Composite Materials 2001;35(10):883— 
904. 


147 


6 Literature 


[84] 


[85] 


[87] 


[88] 


[89] 


148 


Kim YK, White SR. Stress relaxation behavior of 3501-6 epoxy resin 
during cure. Polym. Eng. Sci. 1996;36(23):2852-62. 

Magnus Svanberg J, Anders Holmberg J. Prediction of shape distortions 
Part I. FE-implementation of a path dependent constitutive model. 
Composites Part A: Applied Science and Manufacturing 
2004;35(6):711-21. 

Bogetti TA, Gillespie JW. Process-Induced Stress and Deformation in 
Thick-Section Thermoset Composite Laminates. Journal of Composite 
Materials 1992;26(5):626-60. 

Adolf D, Martin JE. Calculation of Stresses in Crosslinking Polymers. 
Journal of Composite Materials 1996;30(1):13-34. 

Müller V, Kabel M, Andrä H, Böhlke T. Homogenization of linear 
elastic properties of short-fiber reinforced composites — A comparison 
of mean field and voxel-based methods. International Journal of Solids 
and Structures 2015;67-68(8):56—70. 

Müller V, Brylka B, Dillenberger F, Glöckner R, Kolling S, Bohlke T. 
Homogenization of elastic properties of short-fiber reinforced compo- 
sites based on measured microstructure data. Journal of Composite Ma- 
terials 2016;50(3):297-312. 

Voigt W. Ueber die Beziehung zwischen den beiden Elasticitäts- 
constanten isotroper Körper. Ann. Phys. 1889;274(12):573-87. 

Reuss A. Berechnung der Fließgrenze von Mischkristallen auf Grund 
der Plastizitätsbedingung für Einkristalle. Z. angew. Math. Mech. 
1929;9(1):49-58. 

Hashin Z, Shtrikman S. On some variational principles in anisotropic 
and nonhomogeneous elasticity. Journal of the Mechanics and Physics 
of Solids 1962;10(4):335—42. 

Hashin Z, Shtrikman S. A Variational Approach to the Theory of the 
Effective Magnetic Permeability of Multiphase Materials. Journal of 
Applied Physics 1962;33(10):3125-31. 

Hashin Z, Shtrikman S. A variational approach to the theory of the 
elastic behaviour of multiphase materials. Journal of the Mechanics and 
Physics of Solids 1963;11(2):127—40. 

Mori T, Tanaka K. Average stress in matrix and average elastic energy 
of materials with misfitting inclusions. Acta Metallurgica 
1973;21(5):571-4. 


6 Literature 


[96] Tandon GP, Weng GJ. The effect of aspect ratio of inclusions on the 
elastic properties of unidirectionally aligned composites. Polym. Com- 
pos. 1984;5(4):327-33. 

[97] Schapery RA. Thermal Expansion Coefficients of Composite Materials 
Based on Energy Principles. Journal of Composite Materials 
1968;2(3):380-404. 

[98] Rosen BW, Hashin Z. Effective thermal expansion coefficients and 
specific heats of composite materials. International Journal of Engineer- 
ing Science 1970;8(2):157—73. 

[99] Hohberg M. Experimental investigation and process simulation of the 
compression molding process of Sheet Molding Compound (SMC) with 
local reinforcements: Doctoral Thesis, Karlsruher Institut fiir Technolo- 
gie (KIT); 2019. 

[100] Spiess H, Oeckerath A, Landvogt B. Mapping library manual version 
2010.1. Sankt Augustin: Fraunhofer SCAT; 2012. 

[101] Wolf K, Bayrasy P, Brodbeck C, Kalmykov I, Oeckerath A, Wirth N. 
MpCCI: Neutral Interfaces for Multiphysics Simulations. In: Griebel 
M, Schiiller A, Schweitzer MA, editors. Scientific Computing and Al- 
gorithms in Industrial Simulations. Cham: Springer International Pub- 
lishing; 2017, p. 135-151. 

[102] Landervik M, Jergues J. Digimat material model for short fiber rein- 
forced plastics at Volvo Car Corporation. 10th European LS-DYNA 
Conference, Würzburg (Germany) 2015:1-9. 

[103] Reclusado CA, Nagasawa S. Modeling of fiber-reinforced plastics 
taking into account the manufacturing process. In: Papadrakakis M, edi- 
tor. Proceedings of the VII European Congress on Computational 
Methods in Applied Sciences and Engineering (ECCOMAS Congress 
2016). Athens: Institute of Structural Analysis and Antiseismic Re- 
search School of Civil Engineering National Technical University of 
Athens (NTUA) Greece; 2016 - 2016, p. 2220-2235. 

[104] Görthofer J, Meyer N, Pallicity TD, Schöttl L, Trauth A, Schemmann 
M et al. Virtual process chain of sheet molding compound: Develop- 
ment, validation and perspectives. Composites Part B: Engineering 
2019;169(1):133-47. 


149 


6 Literature 


[105] Pipes RB, Hearle JWS, Beaussart AJ, Sastry AM, Okine RK. A Consti- 
tutive Relation for the Viscous Flow of an Oriented Fiber Assembly. 
Journal of Composite Materials 1991;25(9):1204-17. 

[106] Wittemann F, Karger L, Henning F. Theoretical approximation of 
hydrodynamic and fiber-fiber interaction forces for macroscopic simu- 
lations of polymer flow process with fiber orientation tensors. Compo- 
sites Part C: Open Access 2021;132(53):100152. 

[107] Toll S. Packing mechanics of fiber reinforcements. Polym. Eng. Sci. 
1998;38(8):1337-50. 

[108] Toll S, Manson J-AE. Elastic compression of a fibre network. J. Appl. 
Mech. 1995;62:223-6. 

[109] Servais C, Manson J-AE, Toll S. Fiber—fiber interaction in concentrated 
suspensions: Disperse fibers. Journal of Rheology 1999;43(4):991— 
1004. 

[110] Meyer N, Hrymak AN, Karger L. Modeling Short-Range Interactions in 
Concentrated Newtonian Fiber Bundle Suspensions. International Pol- 
ymer Processing 2021;36(3):255-63. 

[111]Pipes RB, Hearle JWS, Beaussart AJ, Okine RK. Influence of Fiber 
Length on the Viscous Flow of an Oriented Fiber Assembly. Journal of 
Composite Materials 1991;25(10):1379-90. 

[112]Pipes RB. Anisotropic Viscosities of an Oriented Fiber Composite with 
a Power-Law Matrix. Journal of Composite Materials 
1992;26(10):1536-52. 

[113] Loredo A, Klöcker H. Generalized inverse of the compliance tensor, 
and behaviour of incompressible anisotropic materials — Application 
to damage. Mechanics Research Communications 1997;24(4):37 1-6. 

[114] Ericksen JL. Transversely isotropic fluids. Kolloid-Zeitschrift 
1960; 173(2):117—22. 

[115] Hinch EJ, Leal LG. Time-dependent shear flows of a suspension of 
particles with weak Brownian rotations. J. Fluid Mech. 
1973;57(04):753. 

[116] Bertöti R, Böhlke T. Flow-induced anisotropic viscosity in short FRPs. 
Mech Adv Mater Mod Process 2017;3(1):751. 

[117] Fan X, Phan-Thien N, Zheng R. A direct simulation of fibre suspen- 
sions. Journal of Non-Newtonian Fluid Mechanics 1998;74(1-3):113- 
35. 


150 


6 Literature 


[118] Lindström SB, Uesaka T. Simulation of the motion of flexible fibers in 
viscous fluid flow. Physics of Fluids 2007;19(11):113307. 

[119] Toll S. Note: On the tube model for fiber suspensions. Journal of 
Rheology 1993;37(1):123-5. 

[120] Férec J, Ausias G, Heuzey MC, Carreau PJ. Modeling fiber interactions 
in semiconcentrated fiber suspensions. Journal of Rheology 
2009;53(1):49-72. 

[121]Toll S, Manson J-AE. Dynamics of a planar concentrated fiber suspen- 
sion with non-hydrodynamic interaction. Journal of Rheology 
1994;38(4):985-97. 

[122] Yamane Y, Kaneda Y, Dio M. Numerical simulation of semi-dilute 
suspensions of rodlike particles in shear flow. Journal of Non- 
Newtonian Fluid Mechanics 1994(54):405—21. 

[123] Bounoua S, Lemaire E, Férec J, Ausias G, Zubarev A, Kuzhir P. Ap- 
parent yield stress in rigid fibre suspensions: the role of attractive col- 
loidal interactions. J. Fluid Mech. 2016;802:611-33. 

[124] Djalili-Moghaddam M, Toll S. A model for short-range interactions in 
fibre suspensions. Journal of Non-Newtonian Fluid Mechanics 
2005;132(1-3):73-83. 

[125] DiBenedetto AT. Prediction of the glass transition temperature of 
polymers: A model based on the principle of corresponding states. J. 
Polym. Sci. B Polym. Phys. 1987;25(9):1949--69. 

[126] Simon SL, Gillham JK. Cure kinetics of a thermosetting liquid dicya- 
nate ester monomer/high-Tg polycyanurate material. J. Appl. Polym. 
Sci. 1993;47(3):461-85. 

[127] Karkanas PI, Partridge IK, Attwood D. Modelling the Cure of a Com- 
mercial Epoxy Resin for Applications in Resin Transfer Moulding. 
Polym. Int. 1996;41(2):183-91. 

[128] Ding A, Li S, Wang J, Ni A, Zu L. A new path-dependent constitutive 
model predicting cure-induced distortions in composite structures. 
Composites Part A: Applied Science and Manufacturing 2017;95:183— 
96. 

[129] Trouton T. On the coefficient of viscous traction and its relation to that 
of viscosity. Proceedings of the Royal Society of London A 
1906;77(519):426-40. 


151 


6 Literature 


[130] Chiba K, Nakamura K. Numerical solution of fiber suspension flow 
through a complex channel. Journal of Non-Newtonian Fluid Mechan- 
ics 1998;78(2-3):167-85. 

[131]sbhpp. Datasheet Porophen GF 9202 L12. [January 11, 2021]; Availa- 
ble from: https://sbhpp- 
static.s3.amazonaws.com/documents/ds_2124_e.pdf. 

[132] Maertens R, Hees A, Schöttl L, Liebig W, Elsner P, Weidenmann KA. 
Fiber shortening during injection molding of glass fiber-reinforced phe- 
nolic molding compounds: fiber length measurement method develop- 
ment and validation. Polymer-Plastics Technology and Materials 
2021:1-14. 

[133] Saalbach H, Maenz T, Englich S, Raschke K, Scheffler T, Wolf S et al. 
Faserverstärkte Duroplaste für die Großserienfertigung im Spritzgießen. 
Ergebnisbericht des BMBF-Verbundprojektes FiberSet 29. Bekanntma- 
chung vom 15.04.2010 "Energieeffizienter Leichtbau" betreut vom 
PTKA Projektträger Karlsruhe im Karlsruher Institut für Technologie: 
[KraussMaffei Technologies GmbH]; 2015. 

[134] Goris S, Back T, Yanev A, Brands D, Drummer D, Osswald TA. A 
novel fiber length measurement technique for discontinuous fiber- 
reinforced composites: A comparative study with existing methods. 
Polym. Compos. 2018;39(11):4058-70. 


152 


Karlsruher Schriftenreihe Fahrzeugsystemtechnik 
FAST Institut für Fahrzeugsystemtechnik 
(ISSN 1869-6058) 


Eine vollständige Übersicht der Bände finden Sie im Verlagsshop 


Band 76 


Band 77 


Band 78 


Band 79 


Band 80 


Band 81 


Band 82 


Band 83 


Kai-Lukas Bauer 

Echtzeit-Strategieplanung für vorausschauendes 
automatisiertes Fahren 

ISBN 978-3-7315-0949-3 


Thomas Schirle 

Systementwurf eines elektromechanischen Fahrwerks 
für Megacitymobilität 

ISBN 978-3-7315-0995-0 


Dominik Dörr 

Simulation of the thermoforming process 

of UD fiber-reinforced thermoplastic tape laminates 
ISBN 978-3-7315-0998-1 


Dominik Robert Naake 

Simulation of damage mechanisms in weave reinforced 
materials based on multiscale modeling 

ISBN 978-3-7315-1005-5 


Martin Hohberg 

Experimental investigation and process simulation of the 
compression molding process of Sheet Molding Compound 
(SMC) with local reinforcements 

ISBN 978-3-7315-1007-9 


Benedikt Fengler 

Manufacturing-constrained multi-objective optimization of 
local patch reinforcements for discontinuous fiber reinforced 
composite parts 

ISBN 978-3-7315-1006-2 


Johannes Masino 

Road Condition Estimation with Data Mining Methods 
using Vehicle Based Sensors 

ISBN 978-3-7315-1004-8 


11. Kolloquium Mobilhydraulik 
10. September 2020, Karlsruhe 
ISBN 978-3-7315-1036-9 


Die Bande sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar. 


Band 84 


Band 85 


Band 86 


Band 87 


Band 88 


Band 89 


Band 90 


Band 91 


Band 92 


Felix Weber 

Beitrag zur Entwicklung von Konstantflusspumpen 
für Frischbeton unter genauerer Betrachtung der 
Dickstoffventile 

ISBN 978-3-7315-1037-6 


8. Fachtagung 

Hybride und energieeffiziente Antriebe für mobile 
Arbeitsmaschinen. 23. Februar 2021, Karlsruhe 
ISBN 978-3-7315-1071-0 


Sebastian Fünfgeld 

Vorausschauende Regelung von Fahrzeugsystemen durch 
stochastische Vorhersage der Fahrzeugdynamik 

ISBN 978-3-7315-1060-4 


Isabelle Charlotte Ays 

Development of a CO2e quantification method and of solutions 
for reducing the greenhouse gas emissions of construction 
machines = Entwicklung einer CO2e Quantifizierungsmethode 
und von Lösungen zur Reduzierung von Treibhausgasemissio- 
nen in Baumaschinen 

ISBN 978-3-7315-1033-8 


Alexander Bernath 

Numerical prediction of curing and process-induced distortion 
of composite structures 

ISBN 978-3-7315-1063-5 


Nils Bulthaupt 

Objektivierung des Schwingungskomforts 
schwerer Nutzfahrzeuge 

ISBN 978-3-7315-1075-8 


Lars Brinkschulte 

Assistenzsysteme zur Reduktion des Schädigungsverhaltens 
von Komponenten einer mobilen Arbeitsmaschine 

ISBN 978-3-7315-1089-5 


Dominik Dörr 

Adaptive Fahrhinweise für ein längsdynamisches 
Fahrerassistenzsystem zur Steigerung der Energieeffizienz 
ISBN 978-3-7315-1090-1 


Jürgen Römer 

Steuerung und Regelung des Lenkradmoments 
durch Nutzung radselektiver Frontantriebe 
ISBN 978-3-7315-1104-5 


Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar. 


Band 93 


Band 94 


Band 95 


Band 96 


Band 97 


Band 98 


Band 99 


Band 100 


Band 101 


Christian Riese 

Werkzeuge und Konzepte für die Untersuchung 
und Entwicklung zukünftiger Kfz-Bremssysteme 
ISBN 978-3-7315-1125-0 


Yaoqun Zhou 
Dynamisches Bremsverhalten des Reifen-Fahrwerk-Systems 
ISBN 978-3-7315-1156-4 


Stefan Haug 

Ganzheitliche Optimierung einer Axialkolbenpumpe durch 
bedarfsangepasste Entlastung tribologischer Kontakte 
ISBN 978-3-7315-1150-2 


Stefan Scheubner 

Stochastic Range Estimation Algorithms for Electric Vehicles 
using Data-Driven Learning Models 

ISBN 978-3-7315-1166-3 


Yusheng Xiang 

Al and loT Meet Mobile Machines: 
Towards a Smart Working Site 
ISBN 978-3-7315-1165-6 


Nils Meyer 

Mesoscale simulation of the mold filling process 
of Sheet Molding Compound 

ISBN 978-3-7315-1173-1 


Christian Timo Poppe 

Process simulation of wet compression moulding 
for continuous fibre-reinforced polymers 

ISBN 978-3-7315-1190-8 


Torben Fischer 

Modellprädiktive Regelung eines innovativen 
Thermomanagement-Systems für batterieelektrische Fahrzeuge 
ISBN 978-3-7315-1199-1 


Florian Wittemann 

Fiber-dependent injection molding simulation 
of discontinuous reinforced polymers 

ISBN 978-3-7315-1217-2 


Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar. 


This work presents novel simulation techniques for injection molding of fiber 
reinforced polymers. Injection molding is one of the most applied processes for 
discontinuous reinforced polymers parts. The process conditions such as filling 
rate, temperature, etc. have a significant influence on the final part properties. 
For an adequate prediction of these properties, a process simulation has to 
depict different aspects, including all process steps, being mold filling, holding, 
in-mold solidification and out-of-mold cooling. To enable a fiber-dependent mold 
filling simulation, an anisotropic forth order viscosity tensor, considering fiber 
orientation, length and volume fraction is used within this work. Novel approach- 
es to approximate the hydrodynamic and fiber-contact forces within the homog- 
enized material are presented. The hydrodynamic forces are used in a novel 
fiber breakage modeling approach to predict the transient fiber length distribu- 
tion during mold filling. Due to the anisotropic fiber flow coupling, the fiber 
breakage directly influences the flow modeling and hence the modeled cavity 
pressure. The fiber breakage and cavity pressure are validated with experimental 
data, showing good agreement. Furthermore, novel approaches for anisotropic 
and cure-dependent material modeling during holding and cooling stage are 
presented, to predict residual stresses and warpage. 


ISSN 1869-6058 
ISBN 978-3-7315-1217-2 


Gedruckt auf FSC-zertifiziertem Papier 


